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FILTERED  MASS  DENSITY  FUNCTION  FOR  SUBGRID  SCALE 
MODELING  OF  TURBULENT  DIFFUSION  FLAMES 


Peyman  Givi1  and  Farhad  A.  Jaberi2 
Department  of  Mechanical  and  Aerospace  Engineering 
State  University  of  New  York  at  Buffalo 
Buffalo,  NY  14260-4400 

Abstract 

The  objective  of  this  research  is  to  improve  and  implement  the  filtered 
density  function  (FDF)  methodology  for  subgrid  scale  (SGS)  modeling,  and 
subsequent  large  eddy  simulation  (LES)  of  turbulent  reacting  flows.  This  work 
was  primarily  concentrated  on  the  followings  three  issues:  (1)  Development  of 
the  velocity  filtered  density  function  (VFDF).  (2)  Development  of  the  velocity- 
scalar  filtered  density  function  (VSFDF).  (3)  Implementation  of  the  scalar  fil¬ 
tered  density  function  (SFDF)  for  LES  of  a  turbulent  methane  jet  flame.  In 
efforts  pertaining  to  (l)-(2),  modeled  transport  equation  for  the  VFDF  and 
VSFDF  were  derived.  These  equation  were  solved  via  a  new  Lagrangian  Monte 
Carlo  scheme.  The  model  predictions  were  compared  with  results  obtained  via 
conventional  LES  closures,  and  were  also  appraised  via  comparison  with  direct 
numerical  simulation  (DNS)  data  of  a  turbulent  mixing  layer.  It  was  shown 
that  the  magnitudes  of  the  subgrid  scale  (SGS)  scalar  fluxes  as  predicted  by 
FDF  were  significantly  larger  than  those  predicted  by  conventional  models,  and 
were  much  closer  to  the  filtered  DNS  results.  In  efforts  pertaining  to  (3),  the 
LES/FMDF  predictions  of  a  methane  jet-flame  were  used  to  construct  the  spa¬ 
tial  and  the  compositional  structure  of  the  flame.  Some  preliminary  results  were 
obtained  and  were  compared  with  experimental  data  obtained  at  the  Combus¬ 
tion  Research  Facility  at  the  Sandia  National  Laboratories.  The  overall  flow 
features  as  predicted  by  the  LES/FMDF  methodology  were  in  good  agreements 
with  experimental  data. 


1  Currently:  Department  of  Mechanical  Engineering,  University  of  Pittsburgh,  Pittsburgh,  PA. 

2 Currently:  Department  of  Mechanical  Engineering,  Michigan  State  University,  East  Lansing, 
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1  Introduction 


The  need  for  advanced  computational  methods  for  analysis  of  turbulent  reacting 
flows  is  obvious  [1].  Large  eddy  simulation  (LES)  is  regarded  as  one  of  the  most 
promising  means  of  simulating  such  flows  [2,3].  Amongst  the  various  LES  strategies, 
the  approach  based  on  the  probability  density  function  (PDF)  has  proven  particularly 
effective  [4-19].  The  formal  means  of  conducting  such  LES  is  by  consideration  of  the 
“filtered  density  function”  (FDF)  which  is  essentially  the  filtered  fine-grained  PDF 
of  the  transport  quantities  [20].  Colucci  et  al.  [8]  developed  a  transport  equation 
for  the  FDF  in  constant  density  turbulent  reacting  flows,  jaberi  et  al.  [9]  extended 
the  methodology  for  LES  of  variable  density  flows  by  consideration  of  the  “filtered 
mass  density  function”  (FMDF)  which  is  the  mass  weighted  FDF.  The  fundamental 
property  of  the  PDF  methods  is  exhibited  by  the  closed  form  nature  of  the  chemical 
source  term  appearing  in  the  transport  equation  governing  the  FDF  (FMDF).  This 
property  is  very  important  as  evidenced  in  several  applications  of  FDF  for  LES  of  a 
variety  of  turbulent  reacting  flows  [8-12,16].  See  Givi  [21]  for  a  recent  review. 


2  Accomplishments 

The  goal  of  this  research  was  to  improve  the  capabilities  of  the  FDF  method  and  to 
implement  it  for  LES  of  chemically  reacting  turbulent  flows.  We  feel  that  we  have 
been  very  successful  in  achieving  the  specific  objectives  of  this  work.  These  objectives 
were: 


1.  Development  and  implementation  of  the  velocity  filtered  density  function 
(VFDF)  for  LES. 

2.  Development  and  implementation  of  the  joint  velocity-scalar  filtered  den¬ 
sity  function  (VSFDF)  for  LES. 

3.  Implementation  of  the  scalar  FDF  (SFDF)  for  LES  of  hydrocarbon  diffu¬ 
sion  flames  and  comparison  with  experimental  data. 

The  efforts  pertaining  to  the  first  two  objectives  are  completed  and  are  fully  described 
in  journal  articles  (included  here  as  Appendix  I  and  Appendix  II).  The  work  pertain¬ 
ing  to  the  third  objective  is  premature  for  publication.  Also,  during  the  work  on 


2 


this  project,  the  PI  was  invited  to  deliver  a  Keynote  Lecture  at  the  Third  AFOSR 
International  Conference  on  DNS  and  LES  (Arlington,  TX,  August  5-9,  2001).  A 
copy  of  the  review  article  on  this  tutorial  lecture  is  given  in  Appendix  III.  The  efforts 
pertaining  to  each  of  these  objectives  are  summarized  in  the  next  three  subsections. 


2.1  Velocity  Filtered  Density  Function  for  Large  Eddy  Sim¬ 
ulation  of  Turbulent  Flows 

In  this  part  of  our  work,  a  methodology  termed  the  “velocity  filtered  density  function” 
(VFDF)  was  developed  and  implemented  for  large  eddy  simulation  (LES)  of  turbu¬ 
lent  flows.  In  this  methodology,  the  effects  of  the  unresolved  subgrid  scales  (SGS) 
were  taken  into  account  by  considering  the  joint  probability  density  function  (PDF) 
of  all  of  the  components  of  the  velocity  vector.  An  exact  transport  equation  was  de¬ 
rived  for  the  VFDF  in  which  the  effects  of  the  SGS  convection  appear  in  closed  form. 
The  unclosed  terms  in  this  transport  equation  were  modeled.  A  system  of  stochastic 
differential  equations  (SDEs)  which  yields  statistically  equivalent  results  to  the  mod¬ 
eled  VFDF  transport  equation  was  constructed.  These  SDEs  were  solved  numerically 
by  a  Lagrangian  Monte  Carlo  procedure  in  which  the  Ito-Gikhman  character  of  the 
SDEs  is  preserved.  The  consistency  of  the  proposed  SDEs  and  the  convergence  of 
the  Monte  Carlo  solution  were  assessed  by  comparison  with  results  obtained  by  an 
Eulerian  LES  procedure  in  which  the  corresponding  transport  equations  for  the  first 
two  SGS  moments  are  solved.  The  VFDF  results  were  compared  with  those  obtained 
via  several  existing  SGS  closures.  These  results  were  also  analyzed  via  a  priori  and 
a  posteriori  comparisons  with  results  obtained  by  direct  numerical  simulation  (DNS) 
of  an  incompressible,  three-dimensional  (3D),  temporally  developing  mixing  layer. 

This  work  is  fully  described  in  Appendix  II. 

2.2  Velocity-Scalar  Filtered  Density  Function  for  Large  Eddy 
Simulation  of  Turbulent  Flows 

In  this  part  of  our  work,  a  methodology  termed  the  “velocity-scalar  filtered  density 
function”  (VSFDF)  was  developed  and  implemented  for  large  eddy  simulation  (LES) 
of  turbulent  flows.  In  this  methodology,  the  effects  of  the  unresolved  subgrid  scales 
(SGS)  were  taken  into  account  by  considering  the  joint  probability  density  function 
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(PDF)  of  the  velocity-scalar  field.  An  exact  transport  equation  was  derived  for  the 
VSFDF  in  which  the  effects  of  the  SGS  convection  and  chemical  reaction  are  closed. 
The  unclosed  terms  in  this  equation  were  modeled.  A  system  of  stochastic  differential 
equations  (SDEs)  which  yields  statistically  equivalent  results  to  the  modeled  VSFDF 
transport  equation  was  constructed.  These  SDEs  were  solved  numerically  by  a  La- 
grangian  Monte  Carlo  procedure  in  which  the  Ito-Gikhman  character  of  the  SDEs  is 
preserved.  The  consistency  of  the  proposed  SDEs  and  the  convergence  of  the  Monte 
Carlo  solution  were  assessed  by  comparison  with  results  obtained  by  a  finite  difference 
LES  procedure  in  which  the  corresponding  transport  equations  for  the  first  two  SGS 
moments  were  solved.  The  VSFDF  results  were  compared  with  those  obtained  via 
other  SGS  closures,  and  all  the  results  were  assessed  via  comparison  with  data  ob¬ 
tained  by  direct  numerical  simulation  (DNS)  of  a  temporally  developing  mixing  layer 
involving  transport  of  a  passive  scalar.  It  was  shown  that  the  values  of  both  the  SGS 
and  the  resolved  components  of  all  second  order  moments  including  the  scalar  fluxes 
are  predicted  well  by  VSFDF.  The  sensitivity  of  the  model’s  (empirical)  constants 
were  assessed  and  it  was  shown  that  the  magnitudes  of  these  constants  are  in  the 
same  range  as  that  typically  employed  in  PDF  methods. 

This  work  is  fully  described  in  Appendix  II. 

2.3  Scalar  Filtered  Density  Function  for  Large  Eddy  Simu¬ 
lation  of  a  Diffusion  Flame 

We  have  started  work  on  the  use  of  the  scalar  filtered  density  function  (SFDF)  method 
for  large  eddy  simulation  (LES)  of  the  piloted  jet  flame  configuration  as  considered 
in  the  experiments  of  the  Combustion  Research  Facility  at  the  Sandia  National  Lab¬ 
oratories  [22-26].  This  flame  has  been  the  subject  of  broad  investigations  by  other 
computational/modeling  methodologies  [25].  In  the  experiments,  three  basic  flames 
are  considered,  identified  by  Flames  D,  E,  and  F.  The  geometrical  configuration  in 
these  flames  is  the  same,  but  the  jet  inlet  velocity  is  varied.  In  Flame  D,  the  fuel  jet 
velocity  is  the  lowest  and  the  flame  is  close  to  equilibrium.  The  jet  velocity  increases 
from  flames  D  to  E  to  F,  with  noticeable  non-equilibrium  effects  in  the  latter  two. 

Some  preliminary  work  was  completed  in  which  the  LES/SFDF  is  used  for  predic¬ 
tions  of  the  Sandia  piloted  jet  flame.  In  these  simulations,  combustion  was  modeled 
via  two  chemistry  models:  (1)  an  equilibrium  model  via  realistic  kinetics,  (2)  a  fi- 
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nite  rate,  single-step  model  for  non-equilibrium  flames.  In  (1),  the  LES/SFDF  was 
employed  in  conjunction  with  equilibrium  methane-oxidation  model.  This  model  was 
enacted  via  “flamelet”  simulations  which  consider  a  laminar  counterflow  (opposed 
jet)  flame  configuration  [27-30].  The  full  methane  oxidation  mechanism  of  the  Gas 
Research  Institute  (GRI)  [31,32]  accounting  for  53  species  and  325  elementary  re¬ 
actions  was  employed.  At  low  strain  rates,  the  flame  is  close  to  equilibrium.  Thus, 
the  thermo-chemical  variables  were  determined  completely  by  the  “mixture  fraction.” 
This  flamelet  library  was  coupled  with  our  LES/SFDF  solver  in  which  transport  of 
the  mixture  fraction  is  considered.  It  is  useful  to  emphasize  here  that  the  PDF  of  the 
mixture  fraction  is  NOT  “assumed”  a  priori  (as  done  in  almost  all  other  flamelet  based 
LES  [33-43]);  rather,  it  is  calculated  explicitly  via  the  SFDF.  In  (2),  the  methane 
oxidation  was  modeled  via  a  finite-rate,  single-step  kinetics  model  [44].  The  system 
of  nonlinear  ordinary  different  equations  representing  reactant  conversion  was  solved 
via  LES/SFDF  for  all  of  the  scalar  variables  (mass  fractions  and  enthalpy).  With 
this,  some  aspects  of  non-equilibrium  chemistry  were  taken  into  account,  albeit  in  an 
idealized  manner. 

With  the  equilibrium  chemistry  model,  the  results  obtained  by  LES/SFDF/flamelet 
were  compared  with  Sandia’s  experimental  data  [22-26].  Up  to  now,  we  have  only 
simulated  Flame  D  and  we  have  observed  good  agreements.  But  more  work  is  required 
to  ensure  the  accuracy  of  our  results.  Simulations  of  flames  E  and  F  have  not  been 
yet  attempted. 

In  closure  of  this  section,  we  would  like  to  state  that  since  our  early  work  on  FDF, 
this  methodology  has  been  used  by  several  other  investigators  ( e.g .  Refs.  [12,16,19]). 
Please  see  Appendix  III  for  a  recent  review  of  current  state  of  progress  in  LES/FDF. 


3  Technology  Transition  and  Transfer 

Our  activities  in  regard  to  this  aspect  of  our  work  are  listed  below: 

1.  Performer:  Dr.  Peyman  Givi,  SUNY  at  Buffalo,  (716)  645-2593  (x.  2320). 

Customer:  Dr.  Paul  E.  DesJardin,  Sandia  National  Laboratories,  Albuquerque,  NM, 
(505)  844-7740. 

Result:  Developed  LES-FDF  methodology  to  be  implemented  by  Sandia  for  com- 
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putational  analysis  of  combustion  processes  in  a  highly  shock  driven  turbulent  flow 
field. 

Application:  As  a  part  of  the  ongoing  LDRD  and  WFO  programs  at  Sandia  in  mod¬ 
eling  and  simulation  of  thermobaric  explosions. 

2.  Performer:  Dr.  Peyman  Givi,  SUNY  at  Buffalo,  (716)  645-2593  (x.  2320). 

Customer:  Dr.  M.S.  Anand,  Rolls-Royce,  Indianapolis,  IN,  (317)  230-2828. 

Result:  Developed  LES-FMDF  methodology  to  be  implemented  by  Rolls-Royce  for 
computational  analysis  of  gas  turbine  combustor  flows. 

Application:  Allows  reliable  prediction  of  the  unsteady  reacting  flowfield  in  low- 
emission  premixed  concepts. 


4  Interaction  with  AFRL 

We  are  very  grateful  to  have  the  opportunity  of  interacting  with  several  of  the  re¬ 
searchers  at  the  Wright  Patterson  AFB.  This  interaction  was  initiated  by  the  PI  (P. 
Givi)  by  contacting  Dr.  Balu  Sekar  at  WFRL/PRTC  at  WPAFB.  As  a  part  of  this 
interaction,  the  PI  visited  the  WPAFB  on  June  12,  2002  and  delivered  a  seminar 
entitled:  “State  of  Progress  in  Large  Eddy  Simulation  of  Turbulent  Combustion.” 
The  seminar  was  followed  by  a  technical  exchange  between  the  PI  and  several  of  the 
researchers  at  various  divisions  of  WPAFB. 


5  Enhancement  of  Technology  and  Education 

With  completion  of  this  research,  we  feel  that  we  have  been  able  to  contribute  to 
maintain  U.S.  leadership  in  a  technology  which  is  of  significant  importance  to  DOD 
and  AFOSR.  This  is  a  very  important  matter,  particularly  when  considering  the 
extent  of  ongoing  concentrated  efforts  in  Europe  and  Asia  in  the  efforts  proposed 
here.  To  appreciate  the  seriousness  of  the  competition,  we  indicate  that  the  number 
of  European  researchers  who  are  currently  involved  in  mathematical  &  computational 
modeling  of  turbulent  reacting  flows  is  significantly  larger  than  that  in  the  U.S.  In  fact, 
considering  only  England,  France  and  Germany,  the  number  of  senior  investigators 
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(not  including  graduate  students)  who  are  working  on  the  specific  research  field  of 
“LES  of  turbulent  combustion”  in  these  three  countries  is  larger  than  that  in  the 
whole  U.S.!  This  indicates  the  wide  recognition  of  the  importance  of  this  research 
field  in  other  countries.  We  feel  that  the  approach  proposed  here  is  very  promising 
in  dealing  with  this  challenging  problem  in  addition  to  enhancing  the  current  state 
of  knowledge  in  turbulent  combustion. 

In  order  to  demonstrate  our  visibility  in  this  research,  here  we  shall  list  all  the  awards 
and  some  of  noticeable  achievements  of  the  personnel  involved  in  this  program. 


5.1  Graduate  Students 

Involvement  of  students  in  research  is  an  issue  which  is  taken  very  seriously  at  our 
University  We  are  committed  to  recruiting  excellent  quality  students  and  involving 
them  in  high  caliber  research.  During  the  past  three  years,  the  following  students 
have  been  supported  under  this  Grant. 

1.  Dr.  Laurent  Y.M.  Gicquel  (Ph.D.  2001).  Currently:  Research  Scientist,  CERFACS, 
Toulouse,  France.  AFOSR  support  is  acknowledged  in  Ph.D.  Dissertation  [45]. 

2.  Mr.  Tomasz  G.  Drozda  (M.S.  2002).  Currently:  Ph.D.  candidate  at  the  University 
of  Pittsburgh.  AFOSR  support  is  acknowledged  in  M.S.  Thesis  [46] 

3.  Mr.  M.  Reza  H.  Sheikhi  (M.S.  2002).  Currently:  Ph.D.  candidate  at  the  University 
of  Pittsburgh. 

5.2  Awards  and  Honors 

1.  Peyman  Givi:  Named  William  Kepler  Whiteford  Chair  Professor  of  Engineering 
at  the  University  of  Pittsburgh  (2002). 

2.  Peyman  Givi:  Promoted  to  University  at  Buffalo  Distinguished  Professor  (2002). 

3.  Peyman  Givi:  Received  Professor  of  the  Year  Award,  Tau  Beta  Pi  Engineering 
Honor  Society,  New  York  Nu  Chapter  (2001-2002). 

4.  Tomasz  G.  Drozda:  Second  Prize  winner  at  the  Graduate  Technical  Paper  Compe¬ 
tition  of  AIAA  Northeast  Regional  Student  Conference,  Tesselater  Polytechnic  Insti¬ 
tute,  Troy,  NY.  Title  of  paper:  “A  Hybrid  Stochastic-Deterministic  Methodology  for 
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Large  Eddy  Simulation  of  Scalar  Mixing  and  Reaction  in  Turbulent  Flows,”  April, 

2002. 

5.  Farhad  A.  Jaberi:  Received  Young  Investigator  Award  of  the  Office  of  Naval 
Research  (2000). 

6.  Farhad  A.  Jaberi:  Received  CAREER  Award  of  the  National  Science  Foundation 

(2000). 

7.  Profile  Featured  in: 

•  Reporter.  13  named  UB  Distinguished  Professor  33  (28),  p.  1,  May  9,  2002. 

•  The  Buffalo  News:  Rolls-Royce  using  UB  technique  to  refine  its  engines,  p.  E4, 
March  19,  2002. 

•  SEAS  News:  LES  results  equal  more  expensive  supercomputer  simulations, 
VIII(l),  p.  6,  Spring  (2002). 

•  Reporter:  New  method  produces  “super”  results.  33(21),  p.  1,  March  14,  2002. 

•  Reporter:  Graduates  bringing  recognition  to  lab,  32(20),  p.  4,  February  15, 

2001. 

•  ASEE  Recruitment  Video,  the  NASA  Langley  Research  Center,  Office  of  Edu¬ 
cation,  Hampton,  VA,  1999-2000. 

5.3  Publications 

In  all  of  the  following  publications,  the  support  from  AFOSR  is  gratefully  acknowl¬ 
edged: 

Invited  Review  Articles: 

1.  P.  Givi,  “A  Review  of  Modern  Developments  in  Large  Eddy  Simulation  of  Tur¬ 
bulent  Reactive  Flows,”  Chapter  in  DNS/LES:  Progress  and  Challenges,  pp.  81-92, 
Editors:  C.  Liu,  L.  Sakell  and  T.  Beutner,  Greyden  Press,  Columbus,  OH,  2001. 
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Refereed  Papers: 


1.  M.R.H.  Sheikhi,  T.G.  Drozda,  P.  Givi,  and  S.B.  Pope,  “Velocity-Scalar  Fil¬ 
tered  Density  Function  for  Large  Eddy  Simulation  of  Turbulent  Flows,”  submitted 
to  Physics  of  Fluids  (2002). 

2.  L.Y.M.  Gicquel,  P.  Givi,  F.A.  Jaberi,  and  S.B.  Pope,  “Velocity  Filtered  Density 
Function  for  Large  Eddy  Simulation  of  Turbulent  Flows,”  Physics  of  Fluids,  14(3), 
1196-1213,  2002. 

3.  L.Y.M.  Gicquel,  P.  Givi,  F.A.  Jaberi,  and  S.B.  Pope,  “Velocity  Filtered  Density 
Function  for  Large  Eddy  Simulation  of  a  Turbulent  Mixing  Layer,  in  DNS/LES: 
Progress  and  Challenges,  pp.  327-334,  Editors:  C.  Liu,  L.  Sakell  and  T.  Beutner, 
Greyden  Press,  Columbus,  OH,  2001. 

Conference  Papers: 

1.  M.R.H.  Sheikhi,  T.G.  Drozda,  P.  Givi,  and  S.B.  Pope,  “Velocity-Scalar  Filtered 
Density  Function  for  Large  Eddy  Simulation  of  Turbulent  Flows,”  submitted  for 
presentation  at  the  55th  Annual  Meeting  of  the  Division  of  Fluid  Dynamics  of  the 
American  Physical  Society,  Dallas,  TX,  November  2002. 

2.  P.  Givi,  L.Y.M.  Gicquel,  F.A.  Jaberi  and  S.B.  Pope,  “PDF  Methods  for  Large 
Eddy  Simulation  of  Turbulent  Reactive  Flows,”  Proceedings  of  IUTAM  Symposium 
on  Turbulent  Mixing  and  Combustion,  pp.  96-97,  Kingston,  Ontario,  Canada,  June 
2-6,  2001. 

3.  L.Y.M.  Gicquel,  P.  Givi,  F.A.  Jaberi  and  S.B.  Pope,  “Velocity  Filtered  Density 
Function  for  Large  Eddy  Simulation  of  Turbulent  Flows,”  Bulletin  of  the  American 
Physical  Society,  45(9),  p.  129,  53rd  Annual  Meeting  of  the  Division  of  Fluid  Dy¬ 
namics  of  the  American  Physical  Society,  Washington,  DC,  Nov.  19-21,  2000. 

4.  F.A.  Jaberi,  S.  James  and  P.  Givi,  “Large  Eddy  Simulations  of  Methane  Jet 
Flames,”  Bulletin  of  the  American  Physical  Society ,  45(9),  p.  72,  53rd  Annual  Meet¬ 
ing  of  the  Division  of  Fluid  Dynamics  of  the  American  Physical  Society,  Washington, 
DC,  Nov.  19-21,  2000. 
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A  methodology  termed  the  “velocity  filtered  density  function”  (VFDF)  is  developed  and 
implemented  for  large  eddy  simulation  (LES)  of  turbulent  flows.  In  this  methodology,  the  effects  of 
the  unresolved  subgrid  scales  (SGS)  are  taken  into  account  by  considering  the  joint  probability 
density  function  of  all  of  the  components  of  the  velocity  vector.  An  exact  transport  equation  is 
derived  for  the  VFDF  in  which  the  effects  of  the  SGS  convection  appear  in  closed  form.  The 
unclosed  terms  in  this  transport  equation  are  modeled.  A  system  of  stochastic  differential  equations 
(SDEs)  which  yields  statistically  equivalent  results  to  the  modeled  VFDF  transport  equation  is 
constructed.  These  SDEs  are  solved  numerically  by  a  Lagrangian  Monte  Carlo  procedure  in  which 
the  Ito-Gikhman  character  of  the  SDEs  is  preserved.  The  consistency  of  the  proposed  SDEs  and  the 
convergence  of  the  Monte  Carlo  solution  are  assessed  by  comparison  with  results  obtained  by  an 
Eulerian  LES  procedure  in  which  the  corresponding  transport  equations  for  the  first  two  SGS 
moments  are  solved.  The  VFDF  results  are  compared  with  those  obtained  via  several  existing  SGS 
closures.  These  results  are  also  analyzed  via  a  priori  and  a  posteriori  comparisons  with  results 
obtained  by  direct  numerical  simulation  of  an  incompressible,  three-dimensional,  temporally 
developing  mixing  layer.  ©  2002  American  Institute  of  Physics.  [DOI:  10.1063/1.1436496] 


L  INTRODUCTION 

The  probability  density  function  (PDF)  approach  has 
proven  useful  for  large  eddy  simulation  (LES)  of  turbulent 
reacting  flows.1715  The  formal  means  of  conducting  such 
LES  is  by  consideration  of  the  “filtered  density  function” 
(FDF)  which  is  essentially  the  filtered  fine-grained  PDF  of 
the  transport  quantities.  In  all  previous  contributions,  the 
FDF  of  the  “scalar”  quantities  is  considered:  Gao  and 
O’Brien,3  Colucci  et  al.',6  Reveillon  and  Vervisch,7  and  Zhou 
and  Pereira12  developed  a  transport  equation  for  the  FDF  in 
constant  density  turbulent  reacting  flows.  Jaberi  etal ,8  ex¬ 
tended  the  methodology  for  LES  of  variable  density  flows  by 
consideration  of  the  “filtered  mass  density  function” 
(FMDF),  which  is  essentially  the  mass  weighted  FDF.  The 
fundamental  property  of  the  PDF  methods  is  exhibited  by  the 
closed  form  nature  of  the  chemical  source  term  appearing  in 
the  transport  equation  governing  the  FDF  (FMDF).  This 
property  is  very  important  as  evidenced  in  several  applica¬ 
tions  of  FDF  for  LES  of  a  variety  of  turbulent  reacting 
flows.6'10,12  However,  since  the  FDF  of  only  the  scalar  quan¬ 
tities  are  considered,  all  of  the  “hydrodynamic”  effects  are 
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modeled.  In  all  previous  LES/FDF  simulations,  these  effects 
have  been  modeled  via  “non-FDF”  methods. 

The  objective  of  the  present  work  is  to  extend  the  FDF 
methodology  to  also  include  the  subgrid  scale  (SGS)  velocity 
vector.  This  is  facilitated  by  consideration  of  the  joint  “ve¬ 
locity  filtered  density  function”  (VFDF).  With  the  definition 
of  the  VFDF,  the  mathematical  framework  for  its  implemen¬ 
tation  in  LES  is  established.  A  transport  equation  is  devel¬ 
oped  for  the  VFDF  in  which  the  effects  of  SGS  convection 
are  shown  to  appear  in  closed  form.  The  unclosed  terms  in 
this  equation  are  modeled  in  a  fashion  similar  to  that  in  the 
Reynolds-averaged  simulation  (RAS)  procedures.  A  La¬ 
grangian  Monte  Carlo  procedure  is  developed  and  imple¬ 
mented  for  numerical  simulation  of  the  modeled  VFDF 
transport  equation.  The  consistency  of  this  procedure  is  as¬ 
sessed  by  comparing  the  first  two  moments  of  the  VFDF 
with  those  obtained  by  the  Eulerian  finite  difference  solu¬ 
tions  of  the  same  moments  transport  equations.  The  results  of 
the  VFDF  simulations  are  compared  with  those  predicted  by 
the  Smagorinsky16  closure,  and  the  “dynamic”  Smagorinsky 
model.17'19  The  VFDF  results  are  also  assessed  via  compari¬ 
sons  with  direct  numerical  simulation  (DNS)  data  of  a  three- 
dimensional  (3D)  temporally  developing  mixing  layer  in  a 
context  similar  to  that  of  Vreman  et  al20 

This  work  deals  with  LES  of  the  velocity  field  in  a  con¬ 
stant  density,  nonreacting  flow.  Consideration  of  the  joint 
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velocity-scalar  FDF  (or  FMDF)  in  variable  density,  chemi¬ 
cally  reacting  flows  will  be  the  subject  of  future  work.  It  is  in 
this  context  that  the  approach  has  its  principal  advantage: 
Convective  transport  (of  momentum  and  species)  is  in  closed 
form. 


=  rdui  ,dpfdXj) + TL(Uj  ydp/dXi),  the  SGS  production  rate  ten- 
sor,  Pij=  -  TL(Ui,uk)d{uj)Lldxk-TL(Uj,uk)d{u>)Lldxk , 
and  the  SGS  dissipation  rate  tensor, 

~2  vtl{  du  i  /  dxk ,  dUj  fdxk). 


II.  FORMULATION 

In  the  mathematical  description  of  incompressible  (unit 
density)  turbulent  flows,  the  primary  transport  variables  are 
the  velocity  vector,  h,(x,/)  (i=  1,2,3),  and  the  pressure, 
p(x,t),  field.  The  equations  which  govern  transport  of  these 
variables  in  space  (jtt)  and  time  (t)  are 


dui  du 

-r-=0, 

dxi  .  dt 


du :  dUiUj 

— ^  +  — —  = 


dp  datj 
dxj  dXj  ' 


For  a  Newtonian  fluid,  the  viscous  stress  tensor  cr^  is  repre¬ 
sented  by 

(du :  du;\ 

=rs!)-  ® 

where  v  is  the  kinematic  viscosity  and  is  assumed  constant. 

Large  eddy  simulation  involves  the  spatial  filtering 
operation 


</(x,r)>i=  J 


f{x',t)g(x’,x)dx', 


where  Q  denotes  the  filter  function,  (/(x,r))L  represents  the 
filtered  value  of  the  transport  variable  /(x,r),  and  /'=/ 
“  {f)i  denotes  the  fluctuations  of  /  from  the  filtered  value. 
We  consider  spatially  and  temporally  invariant  and  localized 
filter  functions,  thus  £7(x',x)^G(x'  -x)  with  the 
properties,21  G(x)  =  G(  — x),  and  /*00G(x) Jx=  1.  Moreover, 
we  only  consider  “positive”  filter  functions24  for  which  all 
the  moments  JZODxmG(x)dx  exist  for  m^O.  The  application 
of  the  filtering  operation  to  the  instantaneous  transport  equa¬ 
tions  yields 

<?(«>  )l  n 


dt  dXi 


d{p)L  [  d(o-ij)L 


drL{Ui,Uj) 


where  TL(ui,Uj)={uiUj)L-(ui)L(uj)i  denotes  the  “gener¬ 
alized  SGS  stresses.”  18  These  stresses  satisfy18 

d  d 

Tl  lTdui,Uj)']  +  —  [(uk)LTL{ui,Uj)] 


-n0-+/y 


In  this  equation,  Tijk=  TL(Ui,Uj,uk) -  v(dldxk)[TL(Uj ,Uj)] 
is  the  SGS  turbulent  transport  tensor  where 

TL(ui,Uj,Uk)  =  (uiUjUk)L-(ui)LTL(Uj,Uk)-(uj)LTL(Ui,Uk) 
~ (uk)LTL{Ui ,«;•)- • 18  The  other  terms  are 
the  SGS  pressure- velocity  scrambling  tensor,  II  (/- 


III.  VELOCITY  FILTERED  DENSITY  FUNCTION  (VFDF) 
A.  Definitions 

The  “velocity  filtered  density  function”  (VFDF),  de¬ 
noted  by  PL ,  is  formally  defined  as 

PL{v,x,t)=  f  e[>,u(x',f)]G(x'-x)dx\ 

J  CO 


e[n,«(x,/)]=  S[v -u(x,t)]= II  #lv j- K,(x,t)], 


where  S  denotes  the  delta  function  and  v  is  the  velocity  state 
vector.  The  term  g[i7,if(x,t)]  is  the  “fine-grained” 
density,11’25’26  and  Eq.  (6)  defines  the  VFDF  as  the  spatially 
filtered  value  of  the  fine-grained  density.  With  the  condition 
of  a  positive  filter  kernel,24  PL  has  all  the  properties  of  the 
PDF.26  For  further  developments,  it  is  useful  to  define  the 
“conditional  filtered  value”  of  the  variable  Q(xj)  by 

(Q(xj)\u(x,t)  =  v)L 

f-ZQ(x'>t)g[v,u(x,it)~\G(x,  —  x)dxr 
■<SW‘' - - •  (7) 

where  {a\P)i  denotes  the  filtered  value  of  a  conditioned  on 
fl.  Equation  (7)  implies 

(i)  for  0(x,O  =  c,  (Q{x,t)\v)L=c, 

(ii)  for  e(x,r)  =  2(u(x,r)),  (Q(xJ)\v)l=Q(v), 

(iii)  Integral  property: 

(Q(x,t))L=  J  ^  ( Q(x,t)\v)LPL(v;x9t)dv ,  (8) 


where  c  is  a  constant,  and  Q(xj)=Q(u(x,t))  denotes  the 
case  where  the  variable  Q  is  completely  described  by  the 
variable  u(xft).  From  these  properties  it  follows  that  the 
filtered  value  of  any  function  of  the  velocity  variable  is  ob¬ 
tained  by  integration  over  the  velocity  space 


(Q(x,0)l=Jo 


Q(v)PL(v;x,t)dv. 


B,  VFDF  transport  equation 

The  exact  transport  equation  for  the  VFDF  is  derived  in 
this  section.  Two  forms  of  this  equation  are  considered  simi¬ 
lar  to  those  previously  developed  in  PDF  methods.27*31  The 
starting  point  is  to  consider  the  time-derivative  of  Eq.  (6) 
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dPL(v;x,t)  [°°  du,(x',t)  <?e[i7,it(x',f)] 


Alternatively,  the  conditional  diffusion  term  in  Eq.  (15) 
can  be  represented  as 


XG(x'-x)dx' . 

This  combined  with  Eq.  (7)  yields 

dPjSy^xj)  d  Iditj  \ 
dt  dvt  \  dt  VJi 


v)  PL{v;x,t) 
I L 


Substituting  Eq.  (1)  into  Eq.  (11)  yields 

<?FL(u;x,f)  _  d  f  /  dutuk  \ 

dt  dvA  \  dxi  V I 


With  the  relation 


d  duvA 


and  decompositions 

VkPL=(Uk)LPL+[vk-(Uk)L}PL > 


<?PL(p;x,f) 
dxk  ’ 


dp  \ 

_  S(p)l  n 

dXi  V) 

daik 

\  „  d(<rik)L 

dxk  V 

)LL~  dxk 

r/r 


d(o-ik)L 

--Wr** 

the  VFDF  transport  equation  becomes 

I,  _/  v  ,p  \Ap)l  dPL 

Dt  dxk[^Vk  AAl)  SXi  dvt 

djo-^L  dPL  |  <?  f  /  dp  \  d{p)A 
dxk  dvi  dVf  \\dXi  VjL  dxt  JPl 

_J_\l!^ik  \  d(aik)L\  | 

dvt  I  \dxkVj  dxk  ]  L  ’ 


where  DIDt=dfdt+(uk)L(d/dxk)  denotes  the  “filtered” 
material  derivative. 

Equation  (15)  is  an  exact  transport  equation  for  the 
VFDF.  The  first  term  on  the  right  hand  side  represents  the 
SGS  convection  of  the  VFDF  in  physical  space  and  is  closed. 
The  second  and  third  terms  (which  are  also  in  closed  form) 
represent  the  convection  in  velocity  space  due  to  the  resolved 
pressure  gradient  and  molecular  diffusion,  respectively.  The 
last  two  terms  are  unclosed  and  denote  convective  effects  in 
the  velocity  space  due  to  SGS  pressure  gradient  and  SGS 
diffusion. 


d  /  darik  \ 

■stKw'H 

d2PL(v\x,t)  d2  [  /  dut  du:  \ 

dxkdxk  dvjdvj  \  dxk  dxk  vJPL(v,x,t)  > 


in  which  the  second  term  on  the  right-hand  side  (rhs)  in¬ 
volves  the  conditional  expected  dissipation.  With  this,  the 
alternate  form  of  the  VFDF  transport  equation  is 

DPL-  9  V t  ,  v  n  ,  d(P)LdPL  ,  d2PL 
Dt  *c4[(°4  +  dxt  dVi  dxkdxk 

.  _£_f / \  _^p}l\d 

dVi  \dXiVlI  dxt  )Pl 


dL  I  dut  dUj  \ 
dVjdllj  \  dxk  dxk  f  Pl  ' 


Equation  (17)  is  another  exact  transport  equation  for  the 
VFDF.  The  first  term  on  the  right-hand  side  represents  the 
SGS  convection  of  the  VFDF  in  physical  space  and  is  closed. 
The  second  term  corresponds  to  the  convection  in  the  veloc¬ 
ity  space  due  to  the  resolved  pressure  gradient.  The  third 
term  represents  molecular  diffusion  of  the  VFDF  in  physical 
space.  The  closure  problem  is  associated  with  the  last  two 
terms. 

C.  Modeled  VFDF  transport  equations 

The  generalized  Langevin  model  (GLM)27,32  is  em¬ 
ployed  for  closure  of  the  VFDF  transport  equation.  Here  we 
introduce  two  modeled  VFDF  equations,  which  are  denoted 
by  “VFDF1”  and  “VFDF2.”  These  are  presented  in  order. 
To  close  Eq.  (17),  VFDF1  is 

d  III dp  \  d{p)A  ' 

d 2  I  dUj  dUj  \ 
dv,dVj  \dxk  dxk  VJPl 

ca~^Gu(vr<uj)^P^+  2C°ed^i 

+  vd(Ui)L  d(Uj)L  d2PL  Jjuik  B2Pl 

dxk  dxk  dVjdVj  dxk  dxkdvt'  1 


To  close  Eq.  (15),  VFDF2  is 

A_\llEL  \  ld(Tik  \  ,  dWikh\  D 

dVi  \\dxivj  dXi  \dxkV)  dxk  ]Pl 


d  1 

3..  [GijiVj- (uy)^)  P/J  +  CqE 


dVidVi' 
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Note  that  these  models  (i.e.,  the  first  two  terms  on  the  right- 
hand  sides  of  Eqs.  (18)  and  (19)  are  the  same,  but  that  they 
model  slightly  different  quantities.  With  this  closure,  the  two 
terms  in  GI;  and  e  jointly  represent  the  SGS  pressure-strain 
and  SGS  dissipation.  These  are  modeled  as11,26 


Gir  ~  ^  2  +  4  c°j  ^  ’  6  =  c°k  /Al  ’  e/*’ 

(20) 

where  a>  is  the  SGS  mixing  frequency,  AL  is  the  filter  width, 
*  =  \rdui^ui)  ts  the  SGS  kinetic  energy,  and  e  =  is  the 
SGS  dissipation  rate. 

With  the  GLM,  the  two  forms  of  the  VFDF  transport 
equation  are 

DPl  .  d  d{p)L  dPL  d1PL 

,  »(“,>£  J’-Pl  ,  e(u,)L  PPl 

dxk  dxk  SVjdVj  dxk  dxkdv{ 


(2,) 


for  VFDF1,  and 
DPl  o 

Dt  dx 


d  ,  4  d(p)L  dPL 


^P l  d 


+  2C°'5^  (22) 

for  VFDF2.  Hereinafter,  Eqs.  (21)  and  (22)  are  referred  to  as 
“VFDF1”  and  “VFDF2,”  respectively.  The  difference  be¬ 
tween  these  two  equations  is  in  the  different  treatment  of  the 
closed  viscous  terms. 


for  VFDF2: 


<%)l  f  d{u^L{U))L 

dt  dXi 

#(p)l  t  d\»j)L  dTL(ut,Uj) 
dXj  dx  idxi  dxt 

d  d 

Jt  [rL(Ui  ,uj)]+  —[(uk)LTL(ui,Uj)] 


^  "f* GfoTii^Uj  ,11^) 4* , U fc) 

d(uj)L  d(ui)L 

-  »i(«i  .*»)  — - TL(uj,uk)  — —  +  C0e Sij . 


It  may  be  seen  that  the  zeroth  and  first  moment  equations  are 
identical  (and  exact);  whereas  the  second  central  moment 
equations  differ  by  the  additional  viscous  term  in  VFDF1 
[Eq.  (24)].  A  comparison  of  these  modeled  equations  with 
Eq.  (5)  shows  that  the  GLM  model  implies 

~  n  i;-  -  ( e  ij  -  fe  <5,y)  =  -  C]  w[  rt(  u < ,  u})  -  f k  <50] , 

C,  =  1  +  |C0.  (27) 

This  is  the  same  as  the  Rotta33  model  as  shown  by  Pope.34 
There  are  two  model  constants  in  the  VFDF  equation.  In 
RAS,  typically34-35  Ce=  1,  and  C0~2.1  (C!  =  4.15).  As 
shown  in  Refs.  27,  34  boundedness  of  the  GLM  coefficients 
C0>0  guarantees  that  the  SGS  stress  is  realizable. 

IV.  EQUIVALENT  STOCHASTIC  SYSTEMS 


D.  Transport  equations  for  moments 

The  zeroth,  first,  and  second  moment  equations  corre¬ 
sponding  to  these  two  formulations  are 

for  VFDF1: 


Huj)l  [  <?(« i)L(Uj)L 
dt  dXj 

rfph  t  &{Ujh  dTL{u,,Uj) 
dx j  V  dXjdXi  dXj 

d  d 

Jt\-'rL{ui,uj)^+—[{uk)LTL{ui,uj)] 


=  ~~fok  Ti(M*  V  ^  [  Ti(  «» ’  J 

d(uj)i 

+  GikTL( Uj , uk)  +  Gjk rL( u i  ,uk)  -  rL(ui , uk) 
d(u 

-TL(uj,uk)—^--  +  C0sSij9  (2‘ 


The  solution  of  the  VFDF  transport  equation  provides  all 
the  statistical  information  pertaining  to  the  velocity  vector. 
The  most  convenient  means  of  solving  this  equation  is  via 
the  Lagrangian  Monte  Carlo  scheme.  The  basis  of  this 
scheme  relies  upon  the  principle  of  equivalent  systems.26,32 
Two  systems  with  different  instantaneous  behaviors  may 
have  identical  statistics  and  satisfy  the  same  PDF  transport 
equation.  In  this  context,  the  general  diffusion  process  is 
considered  via  the  following  system  of  stochastic  differential 
equations  (SDEs):26-31-36'37 

dUi(t)=Mi(mMt)lt)dt+E(X(t)M(t);t)dW'l(t) 

+Fij(X(t),U(t);t)dWj(t),  (28) 

where  X{  and  U{  are  probabilistic  representations  of  x  and  w, 
respectively.  The  coefficients  Df  and  M(  are  the  “drift”  in 
the  phase  space  of  position  and  velocity,  respectively.  The 
terms  B  and  E  are  the  “diffusion”  coefficients  for  physical 
and  velocity  spaces,  respectively;  and  W *  and  Wv(  denote 
independent  Wiener-Levy  processes.38  The  tensor  Ftj  repre¬ 
sents  the  dependency  between  the  velocity  and  physical 
spaces.  This  term  is  needed  to  satisfy  the  Ito  condition  for 


1200  Phys.  Fluids,  Vol.  14,  No.  3,  March  2002 


Gicquel  et  al 


B±  0.  A  comparison  of  the  Fokker-Planck  equation  of  Eq. 
(28)  with  the  modeled  VFDF1  transport  equation,  Eq.  (21) 
yields 


d{p)L  d2{u:)L 

*r+2^+w-w.  *. 


B=yj2v,  E=  'JCqs,  Fij=yjlv- 


d{ui)L 


dxj 


Mu 

(29) 


Therefore,  the  proper  SDEs  which  represent  VFDF1  in  the 
Lagrangian  sense  are 


dXi{t)=Ui(t)dt+  \f2vdWxi(t), 


dUj(t)  = 


S(p)l  ,  „  d2{u,)L 

— - \-2v— — - — 

ox  i  BxkoXk 


+  Gij{Uj{t)-{u])L) 


dt 


+  ^dWvi(t)  +  j2^^dWxj(t).  (30) 

This  stochastic  system  is  the  same  as  that  developed  by 
Dreeben  and  Pope29”31  for  RAS. 

For  VFDF2,  due  to  the  absence  of  diffusion  in  physical 
space  we  must  have  B  =  0.  Therefore,  the  corresponding 
SDEs  are 


dU;(t)  = 


dx{ 


d(v,k)L 

dxk 


+  Gy(^(0-<«y)i) 


dt 


+  4c^dWvi{t). 


(31) 


This  system  is  the  same  as  that  suggested  by  Pope26  and 
Haworth  and  Pope27  for  RAS. 

The  primary  difference  between  the  two  formulations 
VFDF1  and  VFDF2  is  due  to  molecular  effects  in  the  spatial 
diffusion  of  the  VFDF.  This  is  explicitly  included  in  the 
VFDF1  formulation  and  is  also  present  in  the  corresponding 
second  moment  equation.  This  difference  is  expected  to  be 
important  in  flows  where  viscous  effects  are  important;  e.g., 
flow  near  solid  boundaries.29"31  Both  of  these  formulation 
are  considered  in  our  numerical  simulations  as  discussed  be¬ 
low. 


V.  NUMERICAL  SOLUTION  PROCEDURE 

Numerical  solution  of  the  modeled  VFDF  transport 
equation  is  obtained  by  a  Lagrangian  Monte  Carlo  proce¬ 
dure.  The  basis  of  this  procedure  is  the  same  as  that  in 
RAS39"41  and  in  previous  LES/FDF.6’8  But  there  are  some 
subtle  differences  which  are  explained  here.  In  the  Lagrang¬ 
ian  description,  the  VFDF  is  represented  by  an  ensemble  of 
N  statistically  identical  Monte  Carlo  particles.  Each  of  these 
particles  carries  information  pertaining  to  its  velocity 
U(rt)(f)  and  position  X(n)(r),  n=  1,2,..  JV.  This  information  is 
updated  via  temporal  integration  of  Eq.  (28).  The  simplest 
means  of  performing  this  integration  is  via  the  Euler- 
Maruyamma  approximation42 


X”(r,+1)=X?(r,)+D?(r,)Ar+^(^)(A01/2^(^), 

+F^tk)(At)mC](tk),  (32) 

where  D”(tk)=Di(X(n\tk),V^\tky,t),  B<n\tk) 
—  B(X(n)(tk),U(n)(tk);t),...  and  £"(tk),  Cj{tk)  are  indepen¬ 
dent  standardized  Gaussian  random  variables.  This  formula¬ 
tion  preserves  the  Markovian  character  of  the  diffusion 
processes43,44  and  facilitates  affordable  computations. 
Higher-order  numerical  schemes  for  solving  Eq.  (28)  are 
available  42  but  one  must  be  cautious  in  using  them  for  LES.6 
Since  the  diffusion  term  in  Eq.  (28)  strongly  depends  on  the 
stochastic  processes,  the  numerical  scheme  must  be  consis¬ 
tent  with  Ito-Gikhman45,46  calculus.  Equation  (32)  exhibits 
this  property. 

The  statistics  are  evaluated  by  consideration  of  the  en¬ 
semble  of  particles  in  a  “finite  volume”  centered  at  a  spatial 
location.  This  ensemble  provides  “one-time”  statistics.  This 
finite  volume  is  characterized  by  a  cubic  box  of  length  A£. 
This  is  necessary  as,  with  probability  one,  no  particle  will 
coincide  with  the  point  as  considered.32  Here,  a  cubic  box  of 
size  A£  is  used  to  construct  the  ensemble  mean,  variances 
and  covariances  of  the  velocity  vector.  These  values  are  used 
in  the  finite  difference  LES  solver  of  Eq.  (4)  as  described 
below. 

The  SGS  dissipation  rate  and  the  SGS  mixing  frequency 
as  required  in  the  solution  of  the  VFDF  are  evaluated  on  the 
finite  difference  grid  points  and  interpolated  to  the  particle’s 
location.  Ideally,  for  reliable  Eulerian  statistics  and  minimum 
numerical  dispersion,  it  is  desired  to  have  the  size  of  the 
sample  domain  infinitesimally  small  (i.e.,  A£— >0)  and  the 
number  of  particles  within  this  domain  infinitely  large.  That 
is 


PL(v;x,t)  * - VN(v;x,t)=2-  2  S(v-u(n)), 

A^£ — *oo  £  fl  E  A£ 

Af— >0 

(33) 

where  VNe  is  the  Eulerian  PDF  constructed  from  the  particle 
ensemble,  n  e  A£  denotes  the  particles  contained  in  an  en¬ 
semble  box  of  length  A£  centered  at  x;  and  NE  is  the  total 
number  of  particles  within  the  box.  With  a  finite  number  of 
particles,  obviously  a  larger  A£  is  needed.  This  compromise  . 
between  the  statistical  accuracy  and  dispersive  accuracy  im¬ 
plies  that  the  optimum  magnitude  of  A£  cannot,  in  general, 
be  specified  a  priori}1’26  This  does  not  diminish  the  capabil¬ 
ity  of  the  procedure,  but  exemplifies  the  importance  of  the 
parameters  governing  the  statistics. 

To  provide  an  estimate  of  the  proper  A£  size,  a  “point 
estimator”  procedure  is  considered.  With  this  procedure,  the 
mean  values  (the  first  moments  of  the  VFDF)  are  evaluated 
by  ensemble  averaging,  and  spatial  variations  of  these  mean 
values  within  the  box  are  ignored.  With  the  discrete  repre¬ 
sentation  [Eq.  (32)],  the  first  two  moments  in  this  procedure 
are  evaluated  via 
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TABLE  I.  Recapitulation  of  the  VFDF  solution  procedures. 
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Finite  difference 
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(34) 


The  point  estimator  is  obviously  subject  to  both  statistical 
errors  and  dispersive  errors  for  A^O. 

To  determine  the  pressure  field,  the  “mean-field  solver” 
is  based  on  the  “compact  parameter”  finite  difference 
scheme  of  Carpenter.47  This  is  a  variant  of  the  McCormack48 
scheme  in  which  fourth-order  compact  differences  are  used 
to  approximate  the  spatial  derivatives,  and  a  second-order 
symmetric  predictor-corrector  sequence  is  employed  for  time 
discretization.  The  numerical  algorithm  is  a  hyperbolic 
solver  which  considers  a  fully  compressible  flow.  Here,  the 
simulations  are  conducted  with  a  low  Mach  number  (A/ 
=s0.3)  to  minimize' compressibility  effects.  All  the  finite  dif¬ 
ference  operations  are  conducted  on  fixed  and  equally  sized 
grid  points.  The  transfer  of  information  from  these  points  to 
the  location  of  the  Lagrangian  particles  is  conducted  via  in¬ 
terpolation.  A  second-order  (bilinear)  interpolation  scheme  is 
used  for  this  purpose.  The  results  of  previous  work  indicate 
no  significant  improvements  with  the  use  of  higher  order 
interpolation  schemes.6 

The  mean-field  solver  also  determines  the  filtered  veloc¬ 
ity  field.  That  is,  there  is  a  “redundancy”  in  the  determina¬ 
tion  of  the  first  filtered  moments  as  both  the  finite  difference 
and  the  Monte  Carlo  procedures  provides  the  solution  of  this 
field.  This  redundancy  is  actually  very  useful  in  monitoring 
the  accuracy  of  the  simulated  results.  Detailed  discussions 
pertaining  to  this  issue  are  provided  in  Refs.  8,  39-41. 

To  establish  the  consistency  of  the  VFDF  solver,  another 
LES  is  also  conducted  in  which  the  modeled  transport  equa¬ 
tions  for  the  filtered  velocity  and  the  generalized  SGS 
stresses  are  solved  strictly  via  the  finite  difference  scheme. 
These  simulations  are  referred  to  as  LES-FD  and  are  only 


applied  for  the  case  corresponding  to  VFDF2.  That  is,  Eqs. 
(25)  and  (26)  are  considered.  Since  the  SGS  transport  terms 
rL(Ui jUj,uk)  are  unclosed  in  Eq.  (26),  the  values  corre¬ 
sponding  to  these  terms  are  taken  from  the  Monte  Carlo 
solver  and  substituted  in  the  SGS  stress  transport  equations. 
The  attributes  of  all  of  the  scheme  are  summarized  in  Table  I, 
with  further  discussions  in  Refs.  6,  39-41. 

VI.  RESULTS 
A.  Flows  simulated 

Simulations  are  conducted  of  a  two-dimensional  (2D) 
planar  jet,  and  a  3D  temporally  developing  mixing  layer.  The 
jet  flow  simulations  are  conducted  primarily  for  establishing 
the  consistency  of  the  Lagrangian  Monte  Carlo  solver.  For 
this  purpose,  2D  simulations  are  sufficient.  To  analyze  the 
overall  performance  of  the  VFDF  and  to  demonstrate  its  full 
capabilities  and  drawbacks,  3D  simulations  are  required. 

In  the  planar  jet,  a  fluid  issues  from  a  jet  of  width  D  into 
a  co-flowing  stream  with  a  lower  velocity.  The  size  of  the 
domain  in  the  stream  wise  ( x )  and  cross-stream  (y)  directions 
are  O^x^lAD  and  -3.5D^y^3.5D.  The  ratio  of  the  co¬ 
flowing  stream  velocity  to  that  of  the  jet  at  the  inlet  is  kept 
fixed  at  0.5.  A  double-hyperbolic  tangent  profile  is  utilized  to 
assign  the  velocity  distribution  at  the  inlet  plane.  The  forma¬ 
tion  of  the  large  scale  coherent  structures  are  expedited  by 
imposing  low  amplitude  perturbations  at  the  inlet.  In  the  fi¬ 
nite  difference  simulations,  the  characteristic  boundary  con¬ 
dition  procedure  of  Ref.  49  is  used  at  the  inlet,  ffee-shear 
boundary  conditions  are  used  at  the  free-streams  and  the 
pressure  boundary  condition  of  Ref.  50  is  used  at  the  outflow. 

The  temporal  mixing  layer  consists  of  two  parallel 
streams  traveling  in  opposite  directions  with  the  same 
speed.51-53  A  hyperbolic  tangent  profile  is  utilized  to  assign 
the  velocity  distribution  at  the  initial  time.  The  simulations 
are  conducted  for  a  cubic  box,  —  L/2^y^L/2,  0 

where  x,  y,  and  z  denote  the  stream  wise,  the  cross¬ 
stream  and  the  spanwise  directions,  respectively;  and  the 
length,  L  is  specified  such  that  L  =  2Np\u,  where  NP  is  the 
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desired  number  of  successive  vortex  pairings  and  Xu  is  the 
wavelength  of  the  most  unstable  mode  corresponding  to  the 
mean  streamwise  velocity  profile  imposed  at  the  initial 
time.  The  flowfield  is  parameterized  in  a  procedure  some¬ 
what  similar  to  that  by  Vreman  et  al20  The  formation  of  the 
large-scale  structures  are  expedited  through  eigenfunction 
based  initial  perturbations.54,55  This  includes  two- 
dimensional20,52,56  and  three-dimensional52,57  perturbations 
with  a  random  phase  shift  between  the  3D  modes.  This  re¬ 
sults  in  the  formation  of  two  successive  vortex  pairings  and 
strong  three-dimensionality. 

The  flow  variables  are  normalized  with  respect  to  se¬ 
lected  reference  quantities.  In  the  jet  flow,  the  jet  exit  veloc¬ 
ity,  and  the  jet  width  are  the  reference  scales.  In  the  temporal 
mixing  layer,  the  reference  length  is  the  half  initial  vorticity 
thickness,  Lr-Sv(t- 0)/2  {Sv  =  AUl\d{u\)Lldy\mzji,  where 
is  the  Reynolds  averaged  value  of  the  filtered  stream- 
wise  velocity  and  AU  is  the  velocity  difference  across  the 
layer).  The  reference  velocity  is  Ur=AUI2. 

B.  Numerical  specifications 

All  finite  difference  simulations  are  conducted  on 
equally  spaced  grid  points  with  grid  spacings  Ax=Ay 
=  A z  (for  3D)  =  A.  The  resolution  for  LES  of  the  planar  jet 
consists  of  201X101  grid  points.  This  allows  simulations 
with  a  Reynolds  number  Re  =  UrDfv  =  14,000.  The  simula¬ 
tions  of  the  temporal  mixing  layer  are  conducted  on  1933  and 
333  points  for  DNS  and  LES,  respectively.  This  allows  simu¬ 
lations  with  Re=UrLr/v  =  5Q . 

To  filter  the  DNS  data,  a  top-hat  function21  of  the  form 
below  is  used 


G(x'-x)  =  II  G(x? 

i=  1 


g(*;-x,H 


r_L 

Al  1  1  2 


(35) 


0  ’!*;-*/!> 


At' 


in  which  ND  denotes  the  number  of  dimensions,  and  AL 
=  2 A. 58  No  attempt  is  made  to  investigate  the  sensitivity  of 
the  results  to  the  filter  function24  or  the  size  of  the  filter.59 

For  VFDF  simulations  of  the  temporal  mixing  layer,  the 
Monte  Carlo  particles  are  initially  distributed  throughout  the 
computational  region.  For  the  jet  flow,  the  particles  are  sup¬ 
plied  in  the  inlet  region  -  l.75D^y^  1.75D.  As  the  par¬ 
ticles  convect  downstream,  this  zone  distorts  as  it  conforms 
to  the  flow  as  determined  by  the  hydrodynamic  field.  The 
simulation  results  are  monitored  to  ensure  the  particles  fully 
encompass  and  extend  well  beyond  regions  of  nonzero  vor¬ 
ticity  with  an  approximately  uniform  particle  number  den¬ 
sity.  All  simulations  are  performed  with  a  uniform 
“weight”26  of  the  Monte  Carlo  particles.  In  the  temporal 
mixing  layer,  due  to  flow  periodicity  in  the  streamwise  and 
spanwise  directions,  if  the  particle  leaves  the  domain  at  one 
of  these  boundaries  new  particles  are  introduced  at  the  other 
boundary  with  the  same  compositional  values.  In  the  cross¬ 


stream  directions,  the  free- slip  boundary  condition  is  satis¬ 
fied  by  the  mirror-reflection  of  the  particles  leaving  through 
these  boundaries.  In  the  planar  jet,  new  particles  are  intro¬ 
duced  through  the  inlet  boundary  at  a  rate  proportional  to  the 
local  flow  velocity  and  with  a  velocity  makeup  dependent  on 
the  cross-stream  direction  only.  When  the  particles  leave  the 
computational  domain  at  the  outflow,  they  are  no  longer 
tracked.  The  density  of  the  Monte  Carlo  particles  is  deter¬ 
mined  by  the  average  number  of  particles  NE  within  the 
ensemble  domain  of  size  A£XA£(XA£).  The  effects  of 
both  of  these  parameters  are  assessed  to  ensure  the  consis¬ 
tency  and  the  statistical  accuracy  of  the  VFDF  simulations. 

All  results  are  analyzed  both  “instantaneously”  and 
“statistically.”  In  the  former,  the  instantaneous  contours 
(snap-shots)  and  scatter  plots  of  the  variables  of  interest  are 
analyzed.  In  the  latter,  the  “Reynolds-averaged”  statistics 
constructed  from  the  instantaneous  data  are  considered.  In 
the  spatially  developing  flows  this  averaging  procedure  is 
conducted  via  sampling  in  time.  In  the  temporal  mixing 
layer,  the  statistics  are  constructed  by  spatial  averaging  over 
the  x-z  plane  of  statistical  homogeneity.  All  Reynolds  aver¬ 
aged  results  are  denoted  by  an  overbar. 

C.  Consistency  and  convergence  assessments 

The  objective  of  this  section  is  to  demonstrate  the  con¬ 
sistency  of  the  VFDF  formulation  and  the  convergence  of  its 
Monte  Carlo  simulation  procedure.  For  this  purpose,  the  re¬ 
sults  via  VFDF  and  LES-FD  are  compared  against  each 
other.  Since  the  accuracy  of  the  finite  difference  procedure  is 
well-established  (at  least  for  the  first-order  filtered  quanti¬ 
ties),  such  a  comparative  assessment  provides  a  good  means 
of  assessing  the  performance  of  the  Monte  Carlo  solution  of 
the  VFDF.  To  do  so,  the  statistical  results  obtained  from  the 
Monte  Carlo  simulations  of  Eq.  (31)  are  compared  with  the 
finite  difference  solution  of  Eqs.  (25)  and  (26).  Also,  no  at¬ 
tempt  is  made  to  determine  the  appropriate  values  of  the 
model  constants;  the  values  suggested  in  the  literature  are 
adopted34  C0=2.1(C,=4.15)  and  Cs=l. 

In  Fig.  1,  the  instantaneous  contour  plots  of  the  vorticity 
are  shown  as  determined  by  (a)  VFDF2  and  (b)  LES-FD. 
This  figure  provides  a  simple  visual  demonstration  of  the 
consistency  of  the  VFDF2.  Scatter  plots  of  (u)L  vs  (v)L  are 
presented  in  Fig.  2.  The  correlation  and  regression  coeffi¬ 
cients  (denoted,  respectively,  by  p  and  r  on  these  figures)  are 
insensitive  to  A£ .  Figures  3  and  4  show  the  Reynolds  aver-  * 
aged  values  of  the  streamwise  velocity  and  several  compo¬ 
nents  of  the  SGS  stress  tensor  for  several  values  of  A£ ,  with 
N£=40  kept  fixed.  It  is  observed  that  the  first  filtered  mo¬ 
ments  as  obtained  by  VFDF  agree  very  well  with  those  via 
LES-FD  even  for  large  A£  values.  However,  smaller  A£  val¬ 
ues  are  required  for  convergence  of  the  VFDF  predicted  SGS 
stresses  to  those  by  LES-FD.  The  relative  difference  between 
the  L2  norms  of  all  of  the  components  of  the  SGS  tensor  as 
a  function  of  (A£/A)2  is  presented  in  Fig.  5.  Extrapolation 
to  A£=0  shows  that  the  “error”  goes  to  zero  as  A£— >0. 

The  influence  of  NE  on  the  first  two  moments  is  shown 
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(a) 


FIG.  1.  Plot  of  the  vorticity  field  con¬ 
tours,  (a)  VFDF2,  (b)  LES-FD.  A£ 
=  A,  JV£=  40. 


in  Figs.  6  and  7.  It  is  observed  that  NE  does  not  have  a 
significant  influence  on  the  first  moments,  but  does  slightly 
influence  the  second  moments.  In  all  the  cases  considered, 
V£5=40  yields  reliable  predictions,  consistent  with  previous 
consistency  and  convergence  assessments  of  the  scalar 
FDF.6’8  All  the  subsequent  simulations  are  conducted  with 
A£  =  A/2  and  NE=4Q. 

D.  Comparative  assessments  of  the  VFDF 

The  objective  of  this  section  is  to  analyze  some  of  the 
characteristics  of  the  VFDF  via  comparative  assessments 
against  DNS  data.  This  assessment  is  done  via  both  a  priori 
and  a  posteriori  analyses.  In  the  former,  the  DNS  results  are 
used  to  determine  the  range  of  the  empirical  constants  ap¬ 
pearing  in  the  VFDF  sub-closures.  In  the  latter,  the  final  re¬ 
sults  as  predicted  by  the  VFDF  are  directly  compared  with 
those  obtained  by  DNS.  The  procedure  is  similar  to  that  in 
Ref.  20  and  considers  the  3D  temporal  mixing  layer. 

In  addition  to  VFDF,  three  other  LES  are  conducted  with 
(1)  no  SGS  model,  (2)  the  Smagorinsky16,60  SGS  closure, 
and  (3)  the  dynamic  Smagorinsky17-19  model.  In  the  case 
with  no  model,  the  contribution  of  the  SGS  is  completely 
ignored,  i.e.,  rL(ui9Uj)  =  0.  In  this  case,  the  numerical  errors 
amount  to  an  implied  model.  But  as  indicated  in  Ref.  20  this 


case  is  included  to  provide  a  point  of  reference  for  the  other 
closures.  The  Smagorinsky  model  is16’61 

rL(Ui f*  Sij=  -  2  v,S  ij , 


_  ]  (  d(uj)L\ 

"  2  \  dxj  dXi  /’ 


(36) 


*/,=c„a£s. 

Cv=V1  0.172~0.04,  S  =  vV”  and  AL  is  the  characteristic 
length  of  the  filter.  This  model  considers  the  anisotropic  part % 
of  the  SGS  stress  tensor  axj=  rL(ui9Uj)  —  \k 8^ .  The  isotro¬ 
pic  components  are  absorbed  in  the  pressure  field.  The  dy¬ 
namic  version  of  the  Smagorinsky  model  provides  a  means 
of  approximating  Cv  as  suggested  in  Refs.  17-19.  The  pro¬ 
cedure  for  the  implementation  of  this  model  in  the  3D  tem¬ 
poral  mixing  layer  LES  is  described  by  Vreman;20  thus  it  is 
not  repeated  here.  (See  Refs.  11,  23,  62,  and  63  for  recent 
reviews  on  SGS  closure  strategy.) 

In  addition  to  the  resolved  velocity  field,  the  primary 
integral  statistical  quantities  considered  for  comparative  as¬ 
sessments  are 
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FIG.  2.  Scatter  plots  of  the  filtered  velocity  field  as  obtained  via  VFDF2  vs 
LES-FD.  (a),  { u)L ;  (b),  {v)L.  A£=A,  #£=40. 


*/=j 

J 

(37) 

P*-J 

j"  Pfcdx,  with  pk= 

-Tl 

.  d(Ui)L 

» 

Ev=  I 

e„dx,  with  e  = 

v — - - 

M«.>i 

d(Uk)l.\ 

J 

}y  It 

dxk 

\  dxk 

dXi  j’ 

(38) 

B‘=J 

f  min(0,pfc)dx. 

Ef  is  the  kinetic  energy  of  the  resolved  field,  ev  represents 
the  viscous  molecular  dissipation  rate  directly  from  the  fil¬ 
tered  field,  Pk  is  the  production  rate  of  the  SGS  kinetic  en¬ 
ergy  (or  the  rate  of  energy  transfer  from  the  resolved  filtered 
motion  to  the  SGS  motion),  and  Bk  is  the  total 
backscatter. 64-66  The  resolved  molecular  dissipation  rate  is 
always  positive  (by  definition),  but  the  production  rate  of  the 
SGS  kinetic  energy  can  be  locally  negative.  This  backscatter 
is  not  represented  in  the  Smagorinsky  model.  The  dynamic 
model  is  potentially  capable  of  accounting  for  it,  but  at  the 
expense  of  causing  numerical  instabilities.  In  the  implemen¬ 
tation  of  the  dynamic  model  used  here,  backscatter  is 
avoided  by  averaging  the  numerator  and  denominator  of  the 
expression  determining  Cv  (Refs.  19  and  20)  over  the  homo¬ 
geneous  directions.  If  negative  values  are  still  present,  they 
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x 

FIG.  3.  Reynolds  averaged  values  of  the  filtered  streamwise  velocity,  (a) 
Cross-stream  variations  at  x=7,  (b)  streamwise  variation  at  y  =  0  (center- 
line).  Af£=40. 

are  set  equal  to  zero.20,63  The  “resolved”  components  of  the 
Reynolds-averaged  stress  tensor  are  denoted  by  RtJ  where 
^i7=({Mi)t‘“(Mr)z)((M;)jL-(w7)L).  The  “totaP  Reynolds 
stresses  are  denoted  by  r~j  where  =  ( w  — k ,*) ( My  —  uj) . 
These  are  approximated  by  rL(Wj ^y).20’67’68  In 

DNS,  the  total  stresses  are  evaluated  directly  and  the  results 
indicate  that  i?j;+rL(u(-,«;)  does  indeed  approximate  Fjj 
with  a  maximum  error  of  less  than  10%. 

Figure  8  shows  the  distribution  of  the  particle  number 
density  within  the  whole  computational  domain.  Assuring  an 
approximately  uniform  distribution,  the  values  of  the  mo¬ 
ments  within  local  ensembles  are  compared  with  those  of 
filtered  DNS  data.  These  DNS  data  are  transposed  from  the 
original  high  resolution  1933  points  to  the  low  resolution  of 
333  points,  and  then  are  compared  with  LES  results  on  these 
coarse  points. 

The  DNS  data  are  also  used  to  make  a  priori  estimates 
of  the  model  constants.  The  primary  terms  which  require 
closure  are  the  SGS  dissipation  and  the  velocity-pressure 
scrambling  tensors.  The  model  equation  [Eq.  (20)]  involving 
CE  is  in  a  scalar  form.  For  an  estimate  of  Cj  (thus  C0),  we 
consider  the  following  norm  of  the  corresponding  closure 
[Eq.  (27)] 

II  -  ny  -  ( 8  y—  |e  Sij)  || =  Cx  <a\\  tl{  u  ,Uj)  -  \k  <Sy|| ,  (39) 

where  ||  W{j\\  =  yfwJiWy.  To  estimate  the  coefficients,  a  linear 
regression  is  performed  on  all  the  data  points  at  each  com¬ 
putational  time  step.  The  optimized  constants  as  obtained  in 
this  way  are  denoted  by  Ce  and  Cj .  This  procedure  is  also 
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FIG.  4.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  some  of 
the  components  of  the  SGS  stress  tensor  at  x~  7  with  NE~  40.  The  LES-FD 
results  are  obtained  with  A£=0.5A,  NE—AO. 


followed  for  the  Reynolds  averaged  data,  with  the  optimized 
models  obtained  in  this  way  denoted  by  Ce  and  Cj.  The 
temporal  variations  of  these  estimated  values  are  shown  in 
Fig.  9.  The  nonuniformity  of  the  coefficients  indicates  the 
“nonuniversality”  of  the  models.  This  is  expected  as  the  flow 
evolves  from  an  initially  smooth  laminar  state  to  a  strong 
three-dimensional  state  (at  1^40)  before  the  action  of  the 
small  scales  becomes  significant.  The  closures  as  adopted  are 
not  fully  suitable  for  application  in  all  of  these  flow  regions. 
Nevertheless,  Fig.  9  indicates  that  the  values  for  these  coef¬ 
ficients  as  suggested  in  RAS,  i.e.,  Ci~4A5>Ce**l  are  rea¬ 
sonable,  at  least  within  the  turbulent  regime.  The  influences 
of  these  parameters  are  further  investigated  via  a  posteriori 
analysis  of  the  results  as  discussed  below. 

Figures  10  and  11  show  the  contours  of  the  spanwise  and 
the  streamwise  components  of  the  vorticity  field,  respec¬ 
tively,  at  time  t  =  80.  By  this  time,  the  flow  has  gone  through 
several  pairings  and  exhibits  strong  3D  effects.  This  is  evi¬ 
dent  by  the  formation  of  large  scale  spanwise  rollers  with 
presence  of  counter-rotating  streamwise  vortex  pairs  in  all 
the  simulations.  The  results  via  the  no-model  indicate  too 
many  small-scale  structures  which  clearly  are  not  captured 
accurately  on  the  coarse  grid.  The  amount  of  SGS  diffusion 
with  the  Smagorinsky  model  is  very  significant  at  initial 
times.  Due  to  this  dissipative  characteristics  of  the  model,  the 


FIG.  5.  Percentage  of  the  relative  difference  between  the  L2  norms  of  the 
stresses  as  a  function  of  A£/A.  (a)  jc=2.8,  (b)  *=7,  (c)  x=  11.2. 

predicted  results  are  too  smooth  and  only  contain  the  large 
scale  structures.  The  vortical  structures  as  depicted  by  the 
dynamic  Smagorinsky  and  the  VFDF  are  very  similar  and 
predict  the  DNS  results  better  than  the  other  two  models.  The 
results  obtained  by  VFDF1  and  VFDF2  are  virtually  indis¬ 
tinguishable  from  each  other.  This  is  expected,  due  to  the 
lack  of  importance  of  molecular  effects  in  this  free  shear 
flow. 

The  Reynolds  averaged  values  of  the  streamwise  veloc- , 
ity  and  the  temporal  variations  of  the  momentum  thickness 

1  rm  -  - 

<?m(0  4  =  J_t/2(1-<K>t)(1+<«>t)rfy,  (40) 

are  shown  in  Figs.  12  and  13,  respectively.  In  Fig.  12  the 
Reynolds  averaged  values  of  both  filtered  and  unfiltered 
DNS  data  are  considered  and  are  shown  to  be  essentially 
equivalent  Therefore,  the  latter  are  not  shown  in  subsequent 
figures.  The  dissipative  nature  of  the  Smagorinsky  model  at 
initial  times  resulting  in  a  slow  growth  of  the  layer  is  shown. 
Several  values  of  the  model  parameters  (C0,  Ce)  are  consid- 
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FIG.  8.  Particle  number  density  in  VFDF2  simulation  at  f=60.  The  isosur¬ 
face  corresponds  to  NE=A0  set  as  initial  conditions.  C0=2.1,  CB=  1. 


FIG.  6.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  the 
filtered  stream  wise  velocity  at  x=7.  The  LES-FD  results  are  obtained  with 
A£=0.5A,  A^-40. 


LES-FD 
A„=A,Ne=  80 
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ered  in  the  VFDF  simulations.  It  is  observed  that  as  the  mag¬ 
nitude  of  CE  decreases,  the  initial  rate  of  the  layer’s  spread  is 
higher.  With  the  exception  of  the  case  with  Ce  =  0.5  and  the 
Smagorinsky  model,  all  the  other  VFDF  cases,  the  dynamic 
Smagorinsky  and  the  no-model  yield  a  similar  rate  of  layer’s 
growth  at  late  times. 

The  temporal  variations  of  the  resolved  kinetic  energy 
and  all  of  the  terms  defined  in  Eq.  (38)  are  shown  in  Fig.  14. 


t 


FIG.  7.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  some  of 
the  components  of  the  SGS  stress  tensor  at  x  =  7.  The  LES-FD  results  are 
obtained  with  A£~0.5A,  NE=  40. 


FIG.  9.  Time  variation  of  the  model  coefficients  as  obtained  from  a  priori 
analysis  of  the  DNS  data. 
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FIG.  10.  Contour  plots  of  the  spanwise  component  of  the  vorticity  at  z 
=  Q.75L/L, ,  t—  80.  (a)  Filtered  DNS,  (b)  no  model,  (c),  Smagorinsky 
model,  (d)  dynamic  Smagorinsky  model,  (e)  VFDF2,  C0=2.1,  CB=  1,  (f) 
VFDF1,  C0=2.1,  CB=  1. 


The  overall  features  displayed  in  this  figure  are  similar  to 
those  reported  by  Vreman  et  al.20  for  the  no  model,  the  Sma¬ 
gorinsky  model  and  the  dynamic  Smagorinsky  model.  The 
initial  rate  of  decay  of  the  resolved  kinetic  energy  for  the 
Smagorinsky  model  is  the  highest.  This  is  due  to  the  exces¬ 
sive  production  of  the  SGS  kinetic  energy  by  this  model  in 
the  transitional  region,  and  explains  the  reason  for  the  lack  of 
small  scales  in  the  vortical  structures  as  discussed  before.  For 
all  the  other  models  the  initial  rate  of  decrease  of  the  re¬ 
solved  kinetic  energy  is  small  and  increases  as  the  flow  de¬ 
velops.  The  trend  portrayed  by  DNS  results  is  best  captured 
by  the  VFDF  simulations.  For  the  no  model  case  the  only 
means  of  dissipation  of  the  resolved  kinetic  energy  is 
through  molecular  action  and  numerical  dissipation  which 
become  significant  at  later  stages  due  to  presence  of  a  large 
amount  of  small  scales.  In  this  case,  the  amount  of  numerical 
dissipation  is  the  highest.  For  all  the  other  closures,  the  pro¬ 
ductions  rate  of  the  SGS  kinetic  energy  is  larger  than  the 
molecular  dissipation  as  the  flow  develops.  The  dynamic 
Smagorinsky  and  the  no-model  simulations  predict  the  same 
initial  rate  of  decay  for  the  resolved  kinetic  energy.  This  is 
due  to  low  initial  values  of  Pk  predicted  by  the  dynamic 
Smagorinsky  model.  After  t  =  4Q  the  amount  of  Pk  as  pre¬ 
dicted  by  the  dynamic  model  is  more  than  that  of  molecular 
dissipation  by  the  no-model.  Thus  the  rate  of  decay  of  the 
resolved  kinetic  energy  becomes  higher  for  the  dynamic 
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FIG.  11.  Contour  plots  of  the  streamwise  component  of  the  vorticity  vector 
at  jc—  0.25L/Lr,  r-80.  (a)  Filtered  DNS,  (b)  no  model,  (c)  Smagorinsky 
model,  (d)  dynamic  Smagorinsky  model,  (e)  VFDF2,  C0=2.1,  Ce~  1,  (f) 
VFDF1,  C0=2.1,  Cz—  1. 

model  and  is  closer  to  that  obtained  by  DNS. 

With  the  exception  of  the  no-model  case,  all  the  simula¬ 
tions  predict  similar  trends  for  the  molecular  dissipation.  The 
magnitude  of  this  dissipation  as  predicted  by  VFDF  changes 
slightly  with  the  variation  of  the  model  parameter.  The  pro¬ 
duction  rate  of  the  SGS  kinetic  energy  depends  more 
strongly  on  the  model  coefficients;  as  Ce  decreases,  the  peak 


FIG.  12.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  the 
streamwise  velocity  at  f=70. 
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FIG.  13.  Temporal  variations  of  the  momentum  thickness. 


magnitude  of  Pk  is  larger.  The  Smagorinsky  model  does  not 
adequately  predict  Pk ,  and  the  dynamic  model  yields  better 
predictions  at  long  times.  The  overall  trends  are  best  pre¬ 
dicted  by  VFDF.  The  same  is  true  in  capturing  the  backscat- 
ter  phenomenon.  By  design,  the  backscatter  is  identically 
zero  in  the  Smagorinsky  and  the  dynamic  Smagorinsky 
model.  But  VFDF  is  capable  of  capturing  it,  and  its  extent  is 
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FIG.  15.  Cross-stream  variations  of  some  of  the  components  of  ai;  at  t 
=  60. 

controlled  by  the  model  parameters.  In  this  regard  it  is  im¬ 
portant  to  note  that  there  are  no  numerical  instability  prob¬ 
lems  in  the  VFDF  solver  for  negative  Bk  values.  However, 
the  amount  of  predicted  backscatter  is  less  than  that  of 


FIG.  14.  Temporal  variations  of  (a)  total  resolved  ki¬ 
netic  energy,  (b)  SGS  kinetic  energy  production  rate,  (c) 
total  backscatter,  (d)  total  resolved  dissipation. 
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FIG.  16.  Cross-stream  variations  of  some  of  the  components  of  ai}  at  t  FIG.  17.  Cross-stream  variations  of  some  of  the  components  of  Rtj  at  t 
—  80.  =60. 


DNS  and  its  relative  magnitude  is  less  than  those  of  Pk  and 
Ev. 

Several  components  of  the  planar  averaged  values  of  the 
SGS  anisotropy  tensor,  rL(«I  ,m;)- are  presented 
in  Figs.  15  and  16.  Both  the  Smagorinsky  and  the  dynamic 
model  under-predict  the  components  of  this  stress.  The 
VFDF  predictions  are  more  satisfactory.  In  this  regard,  the 
VFDF  is  expected  to  be  more  effective  than  the  other  clo¬ 
sures  for  LES  of  reacting  flows  "since  the  extent  of  SGS  mix¬ 
ing  is  influenced  by  SGS  convection.69,70  “Optimum”  values 
for  Ce  and  C0  cannot  be  suggested  to  predict  all  of  the  com¬ 
ponents  of  this  tensor  at  all  times,  but  it  is  obvious  that  there 
is  too  much  SGS  energy  with  Ce  =  0.5.  _ 

Several  components  of  the  resolved  stress  tensor  are 
shown  in  Figs.  17  and  18.  As  expected,  the  performance  of 
the  Smagorinsky  model  is  not  very  good  as  it  does  not  pre¬ 
dict  the  spread  and  the  peak  value  of  the  resolved  Reynolds 
stresses.  None  of  the  other  models  show  a  distinct  superiority 
in  predicting  the  DNS  results.  The  no-model  and  the  dy¬ 
namic  Smagorinsky  model  predict  large  peak  values  at  the 
middle  of  the  layer.  The  VFDF  predicts  both  the  spread  and 
the  peak  values  reasonably  well.  The  results  for  small  Ce 
values  are  not  shown  since  the  amount  of  energy  in  the  re¬ 


solved  scale  decreases  too  much  in  favor  of  the  increase  of 
the  SGS  stress  (as  shown  in  Figs.  15  and  16).  The  cross¬ 
stream  variations  of  the  total  Reynolds  stress  r12  are  pre¬ 
sented  in  Fig.  19.  The  peak  values  by  the  no-model  simula¬ 
tions  are  again  the  highest.  The  dynamic  model  and  VFDF 
perform  similarly  and  capture  the  DNS  trends  equally  well. 

E.  Comparison  with  previous  investigations 

All  of  the  results  obtained  here  by  DNS,  and  LES  via  the 
Smagorinsky  and  the  dynamic  Smagorinsky  models  agree 
very  well  with  those  of  Vreman  et  al20  The  slight  differences* 
are  due  to  the  nonidentical  flow  initializations,  and  the  dif¬ 
ferent  computational  methodologies  employed  in  the  two 
simulations.  To  compare  with  results  of  other  investigations, 
simulations  are  conducted  of  another  temporally  developing 
mixing  layer  with  Re  =  500  in  a  larger  computational  do¬ 
main,  Lr=  120.  An  initial  forcing  of  the  form  is 

used,  where  A  is  a  uniformly  distributed  random  number 
with  an  amplitude  of  0.05.  Rogers  and  Moser60  perform  DNS 
of  a  high  Re  number  flow  on  512X  210X 192  spectral  points. 
The  results  of  these  simulations  are  in  excellent  agreements 
with  laboratory  data  of  Bell  and  Mehta.71  Here,  LES  is  con 
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FIG.  18.  Cross-stream  variations  of  some  of  the  components  of  R{j  at  t 
=  80. 

ducted  of  this  flow  via  the  dynamic  Smagorinsky  model. 

The  profiles  of  the  mean  streamwise  velocity  and  several 
components  of  the  resolved  stresses  at  t=250  are  presented 
in  Figs.  20  and  21,  respectively.  In  these  figures,  £ 
~y! 4n(0  and  the  symbols  denote  the  experimental  data71  at 
several  streamwise  locations.  The  good  agreement  with  these 
data  also  indicates  good  agreement  with  DNS  results  of  Rog¬ 
ers  and  Moser.72 

F.  Computational  requirements 

The  total  computational  times  associated  with  simula¬ 
tions  of  the  3D  temporal  mixing  layer  are  shown  in  Table  13. 
Expectedly,  the  overhead  associated  with  the  VFDF  simula¬ 
tion  is  extensive  as  compared  to  the  other  models;  neverthe¬ 
less  this  requirement  is  significantly  less  that  of  DNS.  This 
overhead  was  tolerated  in  present  simulations,  but  can  be 
reduced  with  utilization  of  an  optimum  parallel  simulation 
procedure.  This  has  been  discussed  for  use  in  PDF73  and  is 
recommended  for  future  VFDF  simulations. 

VII.  SUMMARY  AND  CONCLUDING  REMARKS 

The  filtered  density  function  (FDF)  methodology1  has 
proven  very  effective  for  large  eddy  simulation  (LES)  of 
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FIG.  19.  Cross-stream  variations  of  r12,  (a)  f  =  60,  (b)  /=80. 

turbulent  reacting  flows.3,6-11  In  all  previous  contributions, 
the  LES/FDF  of  only  the  scalar  quantities  are  considered. 
The  objective  of  the  present  work  is  to  develop  the  FDF 
methodology  for  LES  of  the  velocity  field.  For  this  purpose, 
a  methodology  termed  the  velocity  filtered  density  function 


5 


FIG.  20.  Cross-stream  variation  of  the  Reynolds  averaged  values  of  the 
streamwise  velocity  at  r=250.  Solid  line  denotes  model  predictions  via  the 
dynamic  Smagorinsky  model.  Symbols  denote  experimental  data  of  Bell  and 
Mehta  (Ref.  71). 
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FIG.  21.  Cross-stream  variations  of  the  Reynolds  aver¬ 
aged  values  of  the  streamwise  velocity  at  /=250.  Solid 
lines  denote  model  predictions  via  the  dynamic  Smago- 
rinsky  model.  Symbols  denote  experimental  data  of 
Bell  and  Mehta  (Ref.  71). 


(VFDF)  is  developed.  The  VFDF  is  basically  the  probability 
function  (PDF)  of  the  subgrid  scale  (SGS)  velocity  vector. 
The  exact  transport  equation  governing  the  evolution  of  the 
VFDF  is  derived.  It  is  shown  that  the  effects  of  SGS  convec¬ 
tion  in  this  equation  appears  in  a  closed  form.  The  unclosed 
terms  in  this  transport  equation  are  modeled  via  two  formu¬ 
lations:  VFDF1  and  VFDF2.  The  primary  difference  between 
the  two  models  is  the  inclusion  of  the  molecular  diffusion  in 
the  spatial  transport  of  the  VFDF  in  the  first  formulation.  The 
closure  strategy  in  the  formulation  similar  to  that  in  PDF 
methods  in  Reynolds  averaged  simulation  (RAS) 
procedures.32  In  this  way,  the  VFDF  formulation  is  at  least 
equivalent  to  a  second-order  moment  SGS  closure. 

The  modeled  VFDF  transport  equations  are  solved  nu¬ 
merically  via  a  Lagrangian  Monte  Carlo  scheme  in  which  the 
solutions  of  the  equivalent  stochastic  differential  equations 
(SDEs)  are  obtained.  Two  Monte  Carlo  procedures  are  con¬ 
sidered.  The  schemes  preserve  the  Ito-Gikhman  nature  of 
the  SDEs  and  provide  a  reliable  solution  for  the  VFDF.  The 
consistency  of  the  VFDF  formulation  and  the  convergence  of 
its  Monte  Carlo  solutions  are  assessed.  This  is  done  via  com¬ 
parisons  between  the  results  obtained  by  the  Monte  Carlo 
procedure  and  the  finite  difference  solution  of  the  transport 
equations  of  the  first  two  filtered  moments  of  VFDF  (LES- 
FD).  With  inclusion  of  the  third  moments  from  the  VFDF 
into  the  LES-FD,  the  consistency  and  convergence  of  the 


TABLE  II.  Computer  requirements  for  the  3D  temporal  mixing  layer.  One 
unit  corresponds  to  1657.2  seconds  of  CPU  time  on  the  SGI  origin  2000. 


Resolution 

Ne 

Normalized  CPU  time 

DNS 

193X193X193 

178 

VFDF1 

33X33X33 

40 

33.6 

VFDF2 

33X33X33 

40 

30 

Dynamic  Smagorinsky 

33X33X33 

2.19 

Smagorinsky 

33X33X33 

1.05 

No  model 

33X33X33 

1 

Monte  Carlo  solution  is  demonstrated  by  good  agreements  of 
the  first  two  SGS  moments  with  those  obtained  by  LES-FD. 

The  VFDF  predictions  are  compared  with  those  with 
LES  results  with  no  SGS  model,  with  the  Smagorinsky16 
SGS  closure,  and  with  the  dynamic  Smagorinsky17-19  model. 
All  of  these  results  are  also  compared  with  direct  numerical 
simulation  (DNS)  results  of  a  three-dimensional,  temporally 
developing  mixing  layer  in  a  context  similar  to  that  con¬ 
ducted  by  Vreman  et  al.20  This  comparison  provides  a  means 
of  examining  some  of  the  trends  and  overall  characteristics 
as  predicted  by  LES.  It  is  shown  that  the  VFDF  performs 
well  in  predicting  some  of  the  phenomena  pertaining  to  the 
SGS  transport.  The  magnitude  of  the  SGS  Reynolds  stresses 
as  predicted  by  VFDF  is  larger  than  those  predicted  by  the 
other  SGS  models  and  much  closer  to  the  filtered  DNS  re¬ 
sults.  The  temporal  evolution  of  the  production  rate  of  the 
SGS  kinetic  energy  is  predicted  well  by  VFDF  as  compared 
with  those  via  the  other  closures.  The  VFDF  is  also  capable 
of  accounting  the  SGS  backscatter  without  any  numerical 
instability  problems,  although  the  level  predicted  is  substan¬ 
tially  less  than  that  observed  in  DNS. 

The  results  of  a  priori  assessment  against  DNS  data  in¬ 
dicates  that  the  values  of  the  model  coefficients  as  employed 
in  VFDF  (C0  and  CE)  are  of  the  range  suggested  in  the 
equivalent  models  previously  used  in  RAS.  The  results  of  a 
posteriori  assessments  via  comparison  with  DNS  data  does  * 
not  give  any  compelling  reasons  to  use  values  other  than 
those  suggested  in  RAS,  C0  =  2.1,  C£=  1.  However,  small 
values  of  Cc  are  not  acceptable  as  they  would  yield  too  much 
of  SGS  energy  relative  to  that  within  the  resolved  scales. 

Most  of  the  overall  flow  features,  including  the  mean 
velocity  field  and  the  resolved  and  total  Reynolds  stresses  as 
predicted  by  VFDF  are  similar  to  those  obtained  via  the  dy¬ 
namic  Smagorinsky  model.  This  is  interesting  in  view  of  the 
fact  that  the  model  coefficients  in  VFDF  are  kept  fixed.  It 
may  be  possible  to  improve  the  predictive  capabilities  of  the 
VFDF  by  two  ways:  (1)  Development  of  a  dynamic  proce- 
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dure  to  determine  the  model  coefficients,  and/or  (2)  imple¬ 
mentation  of  higher  order  closures  for  the  generalized  Lange- 
vin  model  parameter  G ^  (see  Ref.  34). 

Work  is  in  progress  towards  developments  of  a  joint 
velocity- sc  alar  FDF  for  LES  of  reacting  flows.  Compared  to 
standard  LES,  this  approach  has  the  advantage  of  treating 
reaction  in  a  closed  form;  and,  compared  to  scalar  PDF6,8  has 
the  advantage  of  treating  convective  transport  (of  momentum 
and  species)  in  closed  form.  These  modeling  advantages 
have  an  associated  computational  penalty.  For  the  cases  con¬ 
sidered  here,  VFDF  is  more  expensive  computationally  than 
the  dynamic  Smagorinsky  model  by  a  factor  of  15.  It  is  ex¬ 
pected  that  VFDF  will  not  be  more  expensive  than  scalar 
FDF,  at  least  for  reacting  flows  with  many  species. 
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Abstract 

A  methodology  termed  the  “velocity-scalar  filtered  density  function”  (VSFDF)  is  developed  and 
implemented  for  large  eddy  simulation  (LES)  of  turbulent  flows.  In  this  methodology,  the  effects 
of  the  unresolved  subgrid  scales  (SGS)  are  taken  into  account  by  considering  the  joint  probability 
density  function  (PDF)  of  the  velocity  and  scalar  fields.  An  exact  transport  equation  is  derived 
for  the  VSFDF  in  which  the  effects  of  the  SGS  convection  and  chemical  reaction  are  closed.  The 
unclosed  terms  in  this  equation  are  modeled  in  a  fashion  similar  to  that  typically  used  in  Reynolds 
averaged  simulation  (RAS)  procedures.  A  system  of  stochastic  differential  equations  (SDEs)  which 
yields  statistically  equivalent  results  to  the  modeled  VSFDF  transport  equation  is  constructed. 
These  SDEs  are  solved  numerically  by  a  Lagrangian  Monte  Carlo  procedure  in  which  the  Ito- 
Gikhman  character  of  the  SDEs  is  preserved.  The  consistency  of  the  proposed  SDEs  and  the 
convergence  of  the  Monte  Carlo  solution  are  assessed  by  comparison  with  results  obtained  by  a 
finite  difference  LES  procedure  in  which  the  corresponding  transport  equations  for  the  first  two  SGS 
moments  are  solved.  The  VSFDF  results  are  compared  with  those  obtained  by  the  Smagorinsky 
model,  and  all  the  results  are  assessed  via  comparison  with  data  obtained  by  direct  numerical 
simulation  (DNS)  of  a  temporally  developing  mixing  layer  involving  transport  of  a  passive  scalar.  It 
is  shown  that  the  values  of  both  the  SGS  and  the  resolved  components  of  all  second  order  moments 
including  the  scalar  fluxes  are  predicted  well  by  VSFDF.  The  sensitivity  of  the  calculations  to  the 
model’s  (empirical)  constants  are  assessed  and  it  is  shown  that  the  magnitudes  of  these  constants 
are  in  the  same  range  as  those  employed  in  PDF  methods. 
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I.  INTRODUCTION 


The  probability  density  function  (PDF)  approach  has  proven  useful  for  large  eddy  sim¬ 
ulation  (LES)  of  turbulent  reacting  flows  [1].  The  formal  means  of  conducting  such  LES 
is  by  considering  the  “filtered  density  function”  (FDF)  [2]  which  is  essentially  the  filtered 
fine-grained  PDF  of  the  transport  quantities.  In  all  previous  contributions,  the  “marginal” 
FDF  of  the  scalars  [3-15],  or  the  marginal  FDF  of  the  velocity  vector  [16]  are  considered; 
see  Givi  [17]  for  a  recent  review. 

The  objective  of  the  present  work  is  to  extend  the  FDF  methodology  to  account  for  the 
“joint”  subgrid  scale  (SGS)  velocity  and  scalar  fields.  This  is  accomplished  by  considering  the 
joint  “velocity-scalar  filtered  density  function”  (VSFDF).  With  the  definition  of  the  VSFDF, 
the  mathematical  framework  for  its  implementation  in  LES  is  established.  A  transport 
equation  is  developed  for  the  VSFDF  in  which  the  effects  of  SGS  convection  and  SGS 
chemical  reaction  (in  a  reacting  flow)  are  closed.  The  unclosed  terms  in  this  equation  are 
modeled  in  a  fashion  similar  to  those  in  the  Reynolds-averaged  simulation  (RAS)  procedures. 
A  Lagrangian  Monte  Carlo  procedure  is  developed  and  implemented  for  numerical  simulation 
of  the  modeled  VSFDF  transport  equation.  The  consistency  of  this  procedure  is  assessed 
by  comparing  the  first  two  moments  of  the  VSFDF  with  those  obtained  by  the  Eulerian 
finite  difference  solutions  of  the  same  moments’  transport  equations.  The  results  of  the 
VSFDF  simulations  are  compared  with  those  predicted  by  the  Smagorinsky  [18]  closure. 
All  the  results  are  assessed  via  comparisons  with  direct  numerical  simulation  (DNS)  data 
of  a  three-dimensional  (3D)  temporally  developing  mixing  layer  involving  transport  of  a 
passive  scalar  variable.  The  sensitivity  of  VSFDF  predictions  to  the  values  of  the  model’s 
(empirical)  constants  is  assessed. 

II.  FORMULATION 

For  the  general  formulation,  we  consider  an  incompressible  (unit  density),  isothermal, 
turbulent  reacting  flow  involving  Ns  species.  The  primary  transport  variables  describing 
such  a  flow  are  the  three  components  of  the  velocity  vector  Uj(x,  t)  (i  =  1,2, 3),  the  pressure 
p(x,  t),  and  the  species’  mass  fractions  (j)a(x,t)  (a  =  1,2  ,...,NS).  The  equations  which 
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govern  the  transport  of  these  variables  in  space  ( Xi )  and  time  (£)  are 


duk 

dxk 

duj  dukut 
dt  dxk 

dfia  .  &'U'k<Pa 
dt  dXk 


=  0 


dp 

dxi 

M 

dxk 


+ 


dojk 

dxk 


+  Sa 


(la) 

(lb) 

(lc) 


where  Sa  =  Sa(<f>(x,t ))  denotes  the  chemical  reaction  term  for  species  a,  and  4>  = 
[4>i,  fa,  ■  ■  ■ ,  <Pns]  denotes  the  scalar  variable  array.  For  an  incompressible,  Newtonian  fluid, 
with  Fick’s  law  of  diffusion,  the  viscous  stress  tensor  <7**,  and  the  scalar  flux  Jk  are  repre¬ 
sented  by 


&ik  — 


duj  duA 

dxk  dxi  J 


4Q  =  -r 


d(pa 

dxk 


(2a) 

(2b) 


where  v  is  the  fluid  kinematic  viscosity  and  T  =  ^  is  the  diffusion  coefficient  of  all  species 
with  Sc  denoting  the  molecular  Schmidt  number.  We  assume  a  constant  value  for  v  =  T;  i.  e. 
Sc  =  1.  In  reactive  flows,  molecular  processes  are  much  more  complicated  than  portrayed 
by  Eq.  (2).  Since  the  molecular  diffusion  is  typically  less  important  than  that  of  SGS,  this 
simple  model  is  adopted  with  justifications  and  caveats  given  by  in  Refs.  [19-21]. 

Large  eddy  simulation  involves  the  spatial  filtering  operation  [1,  22-25] 

/+oo 

f{x',t)G{x',x)dx'  (3) 

•oo 

where  Q(x',x)  denotes  a  filter  function,  and  ( /(x,  t ))  is  the  filtered  value  of  the  transport 
variable  /(x,  t).  We  consider  a  filter  function  that  is  spatially  and  temporally  invariant 
and  localized,  thus:  G(x',x)  =  G(x'  —  x)  with  the  properties  G(x)  >  0,  G(x)dx  =  1. 
Applying  the  filtering  operation  to  Eqs.  (1)  yields 


d  (uk) 
dxk 

djuil  djuk)^{ui) 
dt  dxk 

djda)  d  { Uk )  (<Pa) 
dt  dxk 


=  0 

-  d(P)  |  j d2(^j)  dr(uk,Ui) 
dxi  dxkdxk  dxk 

d2  (</>q)  _  dr^Uk^a)  (s  . 

V  dxkdxk  dxk  “ 


(4a) 

(4b) 

(4c) 
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where  the  second-order  SGS  correlations 


r(a,  b)  =  ( ab )  —  (a)  ( b ) 


(5) 


are  governed  by 

dr^u^Uj)  djuk)  _  db^Uj, Uj) 


dt 


dxk 


2  VT 


......  /  ^(uj)  d(ui ) 

vS^r  ~T{ut'u,y  &T  “ 

dr(uk,Ui,Uj) 


( dui  duj  \  (  dp\  (  dp\ 

\dxfc  ’  dxk )  +  T  V  ’  dxj )  +  T  \  ”  a^  J 


,dxk'  dxk/ 

dT(uh<f)a)  d  (uk)  r(ui,  <fia)  _  d2r(ui,(j)a ) 


<9xi 


+ 


dxk 


=  z/- 


dxkdxk 


2ut 


dp 


+  r  (  —  r(ttj ,  iSq.) 


,  N a(0Q )  ,  x  ^d(ui) 

T{uk,Ui)  —  -  r{uk,(f)Q )-q^~ 

ar(ufc,tfi,0a) 

dxk 


(6a) 


(6b) 


dr{<j)Q,^ ?)  ,  d{uk)r{(pa,(pp )  d2T{(f>a,(pp )  ,  .  xd(<Aa> 

H - O -  =  ^ ^ - T(tzfc,  (pa)—Z - T(Ufc,  0/jJ- - 


dt 


2  i/r 


axfe 
/  d(pa 
\dxk' 


d<b 

dxk 


dxkdxk 

+  x((pa,  Sp)  +  r(<^g,  <?a) 


axfc  axfc 

dx{uk ,  </><*,  </*/j) 
dxk 


(6c) 


In  this  equation,  the  third  order  correlations 

r(o,  6,  c)  =  (a6c)  —  (a)  r(6,  c) 

-  (6)  r(a,  c)  -  (c)  r(a,  6)  -  (a)  (b)  (c)  (7) 

are  unclosed  along  with  the  other  terms  within  square  brackets. 

III.  VELOCITY-SCALAR  FILTERED  DENSITY  FUNCTION  (VSFDF) 

A.  Definitions 

The  “velocity-scalar  filtered  density  function”  (VSFDF),  denoted  by  P,  is  formally  defined 
as  [2] 

/+00 

q ( v ,  ip]  u(x',  f),  <£(x',  t ))  G(x'  -  x)dx'  (8) 

■  OO 

3  Ns 

e  ( V ,  ip]  u(x,  t),  <£(x,  t))  =  5  (^  -  Uj(x,  t))  X  JJ  ^  (^Q  -  ^Q(x,  t ))  (9) 
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where  8  denotes  the  delta  function,  and  v,  ip  are  the  velocity  vector  and  the  scalar  array 
in  the  sample  space.  The  term  g  is  the  “fine-grained”  density  [20,  26],  hence  Eq.  (8) 
defines  VSFDF  as  the  spatially  filtered  value  of  the  fine-grained  density.  With  the  condition 
of  a  positive  filter  kernel  [27],  P  has  all  of  the  properties  of  the  PDF  [20].  For  further 
developments  it  is  useful  to  define  the  “conditional  filtered  value”  of  the  variable  Q(x,  t)  as 


x,  t)  u(x,t) 


(q 

/_r  Q  (Xj>  0  0  K  «(x',  t),  «A(x',  t))  G  (x'  -  x)  dx' 
P(v,1p]X,t) 

(10) 


Equation  (10)  implies  the  following: 


(i)  for  Q(x ,  t)  =  c, 

(ii)  for  Q(x,  t )  =  Q(u(x,  t),  </>(x,  t)) 

(iii)  Integral  properties: 


(<2(x,t)  v,^  =  c  (11a) 

(Q{x,t)  |  v,-0)  =  Q{v,il>)  (lib) 

(Q(x,t))  =  /::  ^Q(x,  t)  |  -u,  i/>^  P(v,  ip]  x,  t)dvdip 


(11c) 


From  Eqs.  (11)  it  follows  that  the  filtered  value  of  any  function  of  the  velocity  and/or  scalar 
variables  is  obtained  by  its  integration  over  the  velocity  and  scalar  sample  spaces 


/+oo 

« 

-oo 


Q(v,  /tp)P(v,  xp]  x,  t)dvdip 


(12) 


B.  VSFDF  Transport  Equations 

To  develop  the  VSFDF  transport  equation,  we  consider  the  time  derivative  of  the  fine¬ 
grained  density  function  (Eq.  (9)) 


dg  _  f  dun  dg  d<pa  dg  \ 
dt  \  dt  dvk  dt  dip a) 


Substituting  Eqs.  (lb-lc),  and  Eqs.  (2a-2b)  into  Eq.  (13)  we  obtain 


dg  dukg 
dt  dxk 


f  dp  v  d2Ui  \  dg_  _  /  d24>a 

\dxi  dxkdxk )  dvi  \  dxkdxk 


+  Sa  (<p) 


dg 

dipa 


(13) 


(14) 
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Integration  of  this  according  to  Eq.  (8),  while  employing  Eq.  (10)  results  in 


dP  dvkP 
dt  dxk 


d  (p)  dP  d 


dxk  dvk  dip a 


+ 


d 

\(/°L 

dvk 

\\dxk 

v , 

7  d2u 

/  it 

i 

dvi  \ 

A  dxkdxk 

[SMP] 

d{p) 


dlpa 


d2(t>a 


dxk 

V,1p\  P 


V- 


v,'tp)P 


(15) 


dxkdxk 

This  is  an  exact  transport  equation  for  the  VSFDF.  It  is  observed  that  the  effects  of  con¬ 
vection  (second  term  on  LHS)  and  chemical  reaction  (the  second  term  on  RHS)  appear  in 
closed  forms.  The  unclosed  terms  denote  convective  effects  in  the  velocity-scalar  sample 
space.  Alternatively,  the  VSFDF  equation  can  be  expressed  as 


dP  dvkP 
dt  dxk 


d2P  d{p)dP  d 


dxkdxk 
d 


+ 


dvk 

d 2 


{{& 


-  2 


dvidvj  [ 

d2 


i ' 


dvidpa 


dxk  dvk 
v,lp 

dui  duj 
dxk  dxk 
dui  d<j), 


dlpa 

d<P) 

dxk 

v,ip  )P 


[SaWP] 


V 


dxk  dxk 


v,i p)P 


_  & 
d<pad(pjs 

This  is  also  an  exact  equation,  but  the  unclosed  terms  are  exhibited  by  the  conditional 
filtered  values  of  the  dissipation  fields  as  shown  by  the  last  three  terms  on  the  RHS. 


v 


d(pg  d(p0 
dxk  dxk 


v,ip)  P 


(16) 
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C.  Modeled  VSFDF  Transport  Equation 


For  closure  of  the  VSFDF  transport  equation,  we  consider  the  general  diffusion  process 
[28],  given  by  the  system  of  stochastic  differential  equations  (SDEs): 

dXt{t)  =  Dx(X+,V+A+\t)dt  +  Bfj(X+,V+,(j>+]t)dWx(t) 

+  FV 7(X+,  U+,  0+;  t)dW?{t)  +  Fx\X+,  U+,  <£+;  t)dWf{t)  (17a) 

dU+(t)  =  Df(X+,U+,<^+;t)dt  +  Bg(X+,U+,^+;t)dWf(t) 

+  F1ujx(X+,  XJ+,<f>+;  t)dWx (t)  +  F^(X+,  U+,  tf>+;  t)dWf{t)  (17b) 
d<j>t{t)  =  Di(x+,v+,(f>+-,t)dt  +  Btj(x+,iJ+A+;t)dwf(t) 

+  F^x(X+,V+^+J)dWx(t)  +  F^(X+,V+A+;t)dWj;(t)  (17c) 

where  X+,U+,il>+  are  probabilistic  representations  of  position,  velocity  vector,  and  scalar 
variables,  respectively.  The  D  terms  denote  drift  in  the  composition  space,  the  B  terms  de¬ 
note  diffusion,  the  F  terms  denote  diffusion  couplings,  and  the  W  terms  denote  the  Wiener- 
Levy  processes  [29,  30].  Following  Haworth  and  Pope  [31],  Dreeben  and  Pope  [32],  Colucci  et 
al.  [7],  and  Gicquel  et  al.  [16]  we  consider  the  generalized  Langevin  model  (GLM)  and  the 
linear  mean  square  estimation  (LMSE)  model  [26] 


dX+  =  U+dt  +  JV^dWx 

dip )  d2  ( u% ) 

dxi  1/2  dxkdxk 


dU +  = 


+  — TT-  +  Gij  {Uf  -  ( Uj )) 


dt 


+  +  y/d*dW? 

oxk 

m  w  -  M + s°w 

+  vvg9-Mdw? 

where  the  variables  u\y  1/2,  * . .  are  all  diffusion  coefficients  (to  be  specified),  and 

=  -uj  +  -C0)  5ij 


dt 


e 

UJ  =  — 

k 


e  =  C, 


kS/2 


k  =  -r  ( uk,uk ) 


(18a) 


(18b) 


(18c) 


(19) 


Here  u>  is  the  SGS  mixing  frequency,  e  is  the  SGS  dissipation  rate,  k  is  the  SGS  kinetic 
energy,  and  A l  is  the  LES  filter  size.  The  parameters  Co,  Cq  and  (\  are  model  constants 
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and  need  to  be  specified.  The  limit  ui  =  uz  =  vsi  =  vs2  —  0  is  the  standard  high  Reynolds 
number  GLM-LMSE  closure  [20]. 

The  Fokker-Planck  equation  [33]  for  f(v,  xp,  x,  t),  the  joint  PDF  of  X+,  U+,  4>+ ,  evolving 
by  the  diffusion  process  as  given  by  Eq.  (18)  is: 


Tt  +  Wk  {Vkf)  ~ 


9W_(t,2.^)s2W 


dxi 


dxkdxk 


1;  [G«  (%-(%»/] 


-  ^  SI \!k+ik  <*■  -  <*«»  n  -  w„  |s°w/l 

.  v\  d2f  - d  (uj)  d2f  - d  (<j>a)  d2/ 

2  dxkdxk  gx.  dxtdvj  vVlVs*  qx.  qx ,Q^a 

,  ^3  d  (ui)  d  (uj)  d2f  1  d2f  ,-—3  (ut)  d  (<f>a)  d2f 
+  —  — - - — — ;; - r  n'-'oe—  ~  +  v^3 VS2 


+ 


2  dxk  dxk  dvidvj  2  dvkdvk 
vs2  d  (<pa)  d  (4>0)  d2f 


dxk  dxk  dvid'ipa 


(20) 


2  dxk  dxk  dipadipp 
The  transport  equations  for  the  filtered  variables  are  obtained  by  integration  of  Eq.  (20) 

according  to  Eq.  (12): 


d(ttfc) 

dxk 

d  ( m )  d  (uk)  (m) 


=  0 


-  +  (I + 1/2  _  IS  “ 


dt 


dxk 


d2  { Ui )  dr(uk)Ui) 


dxk 


d{<f>a)  d{uk)(<f>a)  /  vi\&  (4>a)  (o  dr(uk,<t>a) 

ir+  axt~  +<£»«•)> -—9^7- 


(21a) 


(21b) 


(21c) 


The  transport  equations  for  the  second  order  SGS  moments  are 

dr{ui,Uj)  ,  d{uk)T(ui,Uj)  n  d2r{uuuj)  d  (uj)  8  fa) 

~~di  +  d^k  —  2  dxkdxk  T{Uk’Ui)  dxk  T[Uk’Ui)  dxk 

.  n  -  d(ui)d{uj ) 

+  (yi  —  2y/uiUz  +  Vz) 


dxk  dxk 

+  [GikT{uk,  Uj)  +  Gjkr(uk,  ui)  +  Coe5ij} 


drj^uuUj) 

dxk 


(22a) 


M  ,  9W'(*i,W  "l -SMii,, ft.)  ,  ,S{0a)  ,  ,  ,9W 

at-  + - ail - 7  “  r(“‘’“i)“ST  "  n  *’ w 

+  (yi  -  ~  +  V^2)  g— 


+  [Gifcr(ufc,  (pa)  ~  CnUT(ui,  <pa)]  +  r{uh  SQ)  -  — 


(22b) 
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dr{(j)a,(j)p)  d(uk)T(<f>a,((>0)  V!  d2r{<pa)  fo)  sd  {<t>p}  ,  N  <9  (</>«) 

- 1 - - -  =  — - - - - - T{Uk,(pa)—ZZ - T{Uk,(p0)- - 


dt 


dxk 

+  [yi  -  2sjvxvs2  +  I'Si) 


2  dxkdxk 
d  {(pa)  d  (fo) 


dxk  dxk 

-  [2 CnUT{<f>a,  (p0)]  +  r(0Q,  Sp)  +  t(4>0,  Sa)  - 


dxk 

dr(uk,  0o)  4>p ) 
dxk 


dxk 


(22c) 


A  term-by-term  comparison  of  the  exact  moment  transport  equations  (Eqs.  (4), (6)),  with 
the  modeled  equations  (Eqs.  (21), (22)),  suggests  =  V3  =  i's,  =  vs2  —  2v.  However, 

this  violates  the  realizability  of  the  scalar  field.  A  set  of  coefficients  yielding  a  realizable 
stochastic  model  requires:  z/i  =  i/2  =  1/3  =  2u  and  vs2  =  vs*  =  0.  That  is, 


dX+  = 
dU?  = 


U+dt  +  V^dW? 

d(p)  | 

dxi  dxkdxk 


+  Gy  [u;  -  («,» 


+  \fX?-!X-dwZ  +  VaTedwy 

OXk 

d(f>t  =  -CfiW  (0+  -  (<f>a))  dt 


dt 


(23a) 


(23b) 

(23c) 


The  Fokker-Planck  equation  for  this  system  is 


dt  +  dxk  ^Vkf)  ~  dt!  dli  dvi {Vj  ^  +  dipa  [CnUJ  ^  f] 


+  V 


d2f 


d2f 


,  o  d(uj)  &_L  ,  d{ux)d{u3)  d2f  1  e_ 

dxkdxk  dxi  dxidvj  dxk  dxk  dvidvj  2  0  ~dvkdvk 


(24) 


and  the  corresponding  equations  for  the  moments  are: 


d(uk) 

dxk 

dJVi)  d  ( Uk }  (uj) 

dt  dxk 

dJdoP  djjukjj^} 

dt  dxk 


=  0 

-  d(p)  |  ^.d2(^)  dT(uk,Uj) 

dxi  dxkdxk  dxk 

_  v  —  _  &r(ui'><f>*) 

dxkdxk  dxk 


(25a) 

(25b) 

(25c) 


drj^uUj)  d  (uk)  Tju^Uj)  _  ^rjuiiUj) 
dt  dxk  dxkdxk 


d  {uj}  v  d(ui ) 

r(ufc,«i)— - r(ujk,«j)- 


+  [GifcT(wfc,Uj)  -I-  GjkT(uk,Ui )  +  C0e<5y]  - 


dxk 

dT(uk,Uj,Uj) 

dxk 


dxk 


(26a) 
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(26b) 


dr(ui,  <pa)  d  (uk)  r(ui,  4>a)  d2T{uu<pa)  ,d{<pa)  d  (ui) 

—aT~  + - feT-  =  v~srt ST  ■  T(“‘’U|)-&r  ■  T{Uk'M ^7 


+  [Gifcr(«fc,0a)  -  Cnu)T{uu  <t>a))  +  r(ui ,  Sa)  - 


dr(uk,Ui,(pa) 


dT(<Pon  4>p)  d{uk)T(<pa)<pp)  d2T(<j>a,<t>p)  d((f>p )  d{(f>a ) 

+  — grt — = v  Sxtaxt  -  TiatM-^r  -  T{UkM^r 

+  2^  ~Q~  —  2CctU>T(<f>a,  (ftp)  +  T(<pa,  Sfj)  +  Sa) 

_  drijUk,  fig,  (ftp)  (26c) 

dxk 

which  may  be  compared  to  Eqs.  (4)  and  (6).  Therefore,  the  stochastic  diffusion  process 
described  by  the  SDEs  (23)  implies  the  following  closure  for  the  VSFDF: 

d  \( /  dp  \  djp)\  1  d2  [ / dui  duj  \  ' 
dvk  \\9xfc  V'  )  dxk)  V dvidvj  \dxk  dxk  V'  ) 
d2  [  /  duj  dipa  \  1  d2  [  /  dipa  dj)f}  \  \  ' 

V dvi dipa  \dxkdxk  V'  )  " dip adipp  \dxkdxk\V'  ) 

w  ic  a2/  2i/£M  g/ 

<9xfc  dv^Vj  2  0  dvkdvk  dxk  dxkdvi 

[G«  (»,  -  («,))  /]  +  4L  [Chw  (fc  -  <*,))  /I  (27) 

which  yields  the  closures  at  the  second  order  levels: 


T  ( 0*1’  0*1 )  +  T  (Ui)  0I  )  +  T  (%)  dx)  ~  Gikr{uk' Uj) +  Gjhx(uk,  Ui)  +  C0e5tj 


=  -u;  (  1  + 


r(ui,Uj)  —  ^k5ij  -| eSij  (28a) 


-  2i/r  f— +r  (<Pcn^-)  =GikT(uk,(pa)-Cau}T(ui,<pa) 
\UXk  OXk  J  V  OXi ) 


.^l  +  lCo  +  Cn)rMa)  (28b) 


(28c) 

This  indicates  a  spurious  source  term  in  the  scalar  covariance  equation,  which  is  negligible 
at  high  Reynolds  number  flows. 
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IV.  NUMERICAL  SOLUTION  PROCEDURE 


Numerical  solution  of  the  modeled  VSFDF  transport  equation  is  obtained  by  a  hybrid 
finite  difference-Monte  Carlo  procedure.  The  basis  is  similar  to  those  in  RAS  [34-36]  and 
in  previous  FDF  simulations  [7,  9,  16],  with  some  differences  which  are  described  here. 
For  simulations,  the  FDF  is  represented  by  an  ensemble  of  Np  statistically  identical  Monte 
Carlo  (MC)  particles.  Each  particle  carries  information  pertaining  to  its  position,  X^(i), 
velocity, and  scalar  value,  (p{n\t),  n  —  1, . . . ,  Np.  This  information  is  updated  via 
temporal  integration  of  the  SDEs.  The  simplest  way  of  performing  this  integration  is  via 
Euler-Maruyamma  approximation  [37].  For  example,  for  Eq.  (17a), 

1)  =  Xr(tfc)  +  (Ax(4))nAf+(Rj(4))n(At)^(C/(4))n 

+  {F™{tk))n  (At)1/2  (C»)n  +  (/£♦(**))"  (At)1/2  (c ;/(4))"  (29) 

where  A(4)  =  Di{ytn\tk),U(n\tk),4>(n\tk)-tk),  . . . ,  and  C(4)’s  are  independent  stan¬ 
dardized  Gaussian  random  variables.  This  scheme  preserves  the  Ito  character  of  the  SDEs 

[38]. 

The  computational  domain  is  discretized  on  equally  spaced  finite  difference  grid  points. 
These  points  are  used  for  two  purposes:  (1)  to  identify  the  regions  where  the  statistical 
information  from  the  MC  simulations  are  obtained,  (2)  to  perform  a  set  of  complementary 
LES  primarily  by  the  finite  difference  methodology  for  assessing  the  consistency  and  con¬ 
vergence  of  the  MC  results.  The  LES  procedure  via  the  finite  difference  discretization  is 
referred  to  as  LES-FD  and  will  be  further  discussed  below. 

Statistical  information  is  obtained  by  considering  an  ensemble  of  Ne  computational  par¬ 
ticles  residing  within  an  ensemble  domain  of  characteristic  length  A e  centered  around  each 
of  the  finite-difference  grid  points.  This  is  illustrated  schematically  in  Fig.  1.  For  reliable 
statistics  with  minimal  numerical  dispersion,  it  is  desired  to  minimize  the  size  of  ensem¬ 
ble  domain  and  maximize  the  number  of  the  MC  particles  [20].  In  this  way,  the  ensemble 
statistics  would  tend  to  the  desired  filtered  values: 


«•>* s  k  £  “w 


iV#— >oo 
&E — *0 


(<*> 


te  (“•  *)  =  k  X  (“<n>  -  <“>E)  (iM  -  We)  T  (“. b) 


ti£Ae 


Ne^oo 
Af;— >0 


(30) 
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where  a W  denotes  the  information  carried  by  nth  MC  particle  pertaining  to  transport  vari¬ 
able  a. 

The  LES-FD  solver  is  based  on  the  compact  parameter  finite  difference  scheme  [39,  40]. 
This  is  a  variant  of  the  MacCormack  scheme  in  which  fourth-order  compact  differenc¬ 
ing  schemes  are  used  to  approximate  the  spatial  derivatives,  and  second-order  symmetric 
predictor-corrector  sequence  is  employed  for  time  discretization.  All  of  the  finite  difference 
operations  are  conducted  on  fixed  grid  points.  The  transfer  of  information  from  the  grid 
points  to  the  MC  particles  is  accomplished  via  a  second-order  interpolation.  The  transfer  of 
information  from  the  particles  to  the  grid  points  is  accomplished  via  ensemble  averaging  as 
described  above. 

The  LES-FD  procedure  determines  the  pressure  field  which  is  used  in  the  MC  solver. 
The  LES-FD  also  determines  the  filtered  velocity  and  scalar  fields.  That  is,  there  is  a 
“redundancy”  in  the  determination  of  the  first  filtered  moments  as  both  the  LES-FD  and 
the  MC  procedures  provides  the  solution  of  this  field.  This  redundancy  is  actually  very  useful 
in  monitoring  the  accuracy  of  the  simulated  results  as  shown  in  previous  work  [9,  16,  34-36]. 
To  establish  consistency  and  convergence  of  the  MC  solver,  the  modeled  transport  equations 
for  the  generalized  second  order  SGS  moments  (Eq.  (26)  are  also  solved  via  LES-FD.  In  doing 
so,  the  unclosed  third  order  correlations  are  taken  from  the  MC  solver.  The  comparison  of 
the  first  and  second  order  moments  as  obtained  by  LES-FD  with  those  obtained  by  the  MC 
solver  is  useful  to  establish  the  accuracy  of  the  MC  solver.  These  simulations  are  referred 
to  as  VSFDF-C.  Attributes  of  all  the  simulation  procedures  are  summarized  in  Table  1.  In 
this  table  and  hereinafter,  VSFDF  simulations  refer  to  the  hybrid  MC/LES-FD  procedure 
in  which  the  LES-FD  is  used  for  only  the  first  order  filtered  variables.  In  VSFDF-C,  the 
LES-FD  procedure  is  used  for  both  first  and  second  order  filtered  values.  Further  discussions 
about  the  simulation  methods  are  available  in  Refs.  [7,  16,  34-36]. 

V.  RESULTS 

A.  Flows  Simulated 

Simulations  are  conducted  of  a  two-dimensional  (2D)  and  a  3D  temporally  developing 
mixing  layer  involving  transport  of  a  passive  scalar  variable.  The  2D  simulations  are  per- 
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formed  to  establish  and  demonstrate  the  consistency  of  the  MC  solver.  The  3D  simulations 
are  used  to  assess  the  overall  predictive  capabilities  of  the  VSFDF  methodology.  These 
predictions  are  compared  with  data  obtained  by  direct  numerical  simulation  (DNS)  of  the 
same  layer. 

The  temporal  mixing  layer  consists  of  two  parallel  streams  travelling  in  opposite  direc¬ 
tions  with  the  same  speed  [41-43].  In  the  representation  below,  x,  y  (and  z)  denote  the 
streamwise,  the  cross-stream,  (and  the  span- wise)  directions  (in  3D),  respectively.  The  ve¬ 
locity  components  along  these  directions  are  denoted  by  u,  v,  (and  w)  in  the  x,  y,  (and 
z)  directions,  respectively.  Both  the  filtered  streamwise  velocity  and  the  scalar  fields  are 
initialized  with  a  hyperbolic  tangent  profiles  with  ( u )  =  0.5,  (4>)  =  1  on  the  top  stream  and 
( u )  =  —0.5,  ((f))  =  0  on  the  bottom  stream.  The  length  L  is  specified  such  that  L  =  2Np\u, 
where  Np  is  the  desired  number  of  successive  vortex  pairings  and  Au  is  the  wavelength  of 
the  most  unstable  mode  corresponding  to  the  mean  streamwise  velocity  profile  imposed  at 
the  initial  time.  The  flow  variables  are  normalized  with  respect  to  the  half  initial  vorticity 
thickness,  Lr  =  ,  (5V  =  ■  A.^ — ,  where  (u)L  is  the  Reynolds  averaged  value  of  the 

filtered  streamwise  velocity  and  A U  is  the  velocity  difference  across  the  layer).  The  reference 
velocity  is  Ur  =  AU/2. 

All  2D  simulations  are  conducted  for  0  <  x  <  L,  and  —  y  <  y  <  The  formation 
of  large  scale  structures  is  facilitated  by  introducing  small  harmonic,  phase-shifted,  dis¬ 
turbances  containing  sub-harmonics  of  the  most  unstable  mode  into  the  stream-wise  and 
cross-stream  velocity  profiles.  For  Np  =  1,  this  results  in  formation  of  two  large  vortices 
and  one  subsequent  pairing  of  these  vortices.  The  3D  simulations  are  conducted  for  a  cubic 
box,  0<a;<L,  ^  <  y  <  j,  (0  <  z  <  L).  The  3D  field  is  parameterized  in  a  procedure 
somewhat  similar  to  that  by  Vreman  et  al.  [44].  The  formation  of  the  large  scale  structures 
are  expedited  through  eigenfunction  based  initial  perturbations  [45,  46].  This  includes  two- 
dimensional  [42,  44,  47]  and  three-dimensional  [42,  48]  perturbations  with  a  random  phase 
shift  between  the  3 D  modes.  This  results  in  the  formation  of  two  successive  vortex  pairings 
and  strong  three-dimensionality. 
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B.  Numerical  Specifications 


Simulations  are  conducted  on  equally-spaced  grid  points  with  grid  spacings  Ax  =  Ay  = 
Az  (for  3D)  =  A.  All  2D  simulations  are  performed  on  32  x  41  grid  points.  The  3D 
simulations  are  conducted  on  1933  and  33s  points  for  DNS  and  LES,  respectively.  The 
Reynolds  number  is  Re  =  =  50.  To  filter  the  DNS  data,  a  top-hat  function  of  the  form 

below  is  used 


No  attempt  is  made  to  investigate  the  sensitivity  of  the  results  to  the  filter  function  [27]  or 
the  size  of  the  filter  [49]. 

The  MC  particles  are  initially  distributed  throughout  the  computational  region.  All  sim¬ 
ulations  are  performed  with  a  uniform  “weight”  [20]  of  the  particles.  Due  to  flow  periodicity 
in  the  streamwise  (and  span- wise  in  3D)  direction(s),  if  the  particle  leaves  the  domain  at 
one  of  these  boundaries  new  particles  are  introduced  at  the  other  boundary  with  the  same 
velocity  and  compositional  values.  In  the  cross-stream  directions,  the  free-slip  boundary 
condition  is  satisfied  by  the  mirror-reflection  of  the  particles  leaving  through  these  bound¬ 
aries.  The  density  of  the  MC  particles  is  determined  by  the  average  number  of  particles  Ne 
within  the  ensemble  domain  of  size  A#  x  A#  (x  Ag).  The  effects  of  both  of  these  parameters 
are  assessed  to  ensure  the  consistency  and  the  statistical  accuracy  of  the  VSFDF  simula¬ 
tions.  All  results  are  analyzed  both  “instantaneously”  and  “statistically.”  In  the  former, 
the  instantaneous  contours  (snap-shots)  and  scatter  plots  of  the  variables  of  interest  are  an¬ 
alyzed.  In  the  latter,  the  “Reynolds-averaged”  statistics  constructed  from  the  instantaneous 
data  are  considered.  These  are  constructed  by  spatial  averaging  over  x  (and  z  in  3D).  All 
Reynolds  averaged  results  are  denoted  by  an  overbar. 


C.  Consistency  and  Convergence  Assessments 

The  objective  of  this  section  is  to  demonstrate  the  consistency  of  the  VSFDF  formulation 
and  the  convergence  of  its  MC  simulation  procedure.  For  this  purpose,  the  results  via  MC 


and  LES-FD  are  compared  against  each  other  in  VSFDF-C  simulations.  Since  the  accuracy 
of  the  FD  procedure  is  well-established  (at  least  for  the  first  order  filtered  quantities),  such 
a  comparative  assessment  provides  a  good  means  of  assessing  the  performance  of  the  MC 
solution.  No  attempt  is  made  to  determine  the  appropriate  values  of  the  model  constants; 
the  values  suggested  in  the  literature  are  adopted  [50]  Co  =  2.1,  C£  =  1  and  Co,  =  1.  The 
influence  of  these  parameters  are  assessed  in  Section  (VD). 

The  uniformity  of  the  MC  particles  is  checked  by  monitoring  their  distributions  at  all 
times,  as  the  particle  number  density  must  be  proportional  to  fluid  density.  The  Reynolds 
averaged  density  field  as  obtained  by  both  LES-FD  and  by  MC  are  shown  in  Fig.  2.  Close  to 
unity  values  for  the  density  at  all  times  is  the  first  measure  of  the  accuracy  of  simulations. 
Figures  3  and  4  show  the  instantaneous  contour  plots  of  the  filtered  scalar  and  vorticity 
fields  at  several  times.  These  figures  provide  a  visual  demonstration  of  the  consistency  of 
the  VSFDF.  This  consistency  is  observed  for  all  first  order  moments  without  any  statistical 
variability.  Also,  all  of  these  moments  show  very  little  dependence  on  the  values  of  A#  and 
Ne  consistent  with  previous  FDF  simulations  [7,  9,  16].  In  the  presentation  below  we  only 
focus  on  second  order  moments.  Specifically,  the  scalar-velocity  correlations  are  shown  since 
all  other  second  order  SGS  moments  behave  similarly. 

Figure  5-6  show  the  statistical  variability  of  the  results  for  simulations  with  Ne  =  40. 
It  is  observed  that  these  moments  exhibit  spreads  with  variances  decreasing  as  the  size  of 
the  ensemble  domain  is  reduced.  Figures  7-10  show  the  sensitivity  to  Ne  and  A e-  All 
these  results  clearly  display  convergence  suggested  by  Eq.  (30).  As  the  ensemble  domain 
size  decreases,  the  VSFDF  results  converge  to  those  of  LES-FD.  Ideally,  the  LES-FD  results 
should  become  independent  of  the  MC  results,  as  the  latter  become  more  reliable,  i.e.  when 
Ne  —>  oo,  A e  — »  0).  It  is  observed  that  best  match  is  achieved  with  A e  <  A/2  and 
Ne  >  40.  This  conclusion  is  consistent  with  previous  assessment  studies  on  the  scalar 
FDF  [7,  9],  and  the  velocity  FDF  [16].  All  the  subsequent  simulations  are  conducted  with 
A e  =  A/2  and  Ne  =  40. 

D.  Comparative  Assessments  of  the  VSFDF 

The  objective  of  this  section  is  to  analyze  some  of  the  characteristics  of  the  VSFDF  via 
comparative  assessments  against  DNS  data.  In  addition,  comparisons  are  also  made  with 
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LES  via  the  “conventional”  Smagorinsky  [18,  51]  model: 


T£»(^ij  Uj)  g  k  & ij  2  t'f  Sij , 

1  /  d(ui)L  d(Uj)L\ 
ij  ~  2  V  dxj  dxx  ) 


vt  =  Cu  A|  5,  rt 


Set 


(32) 


C„  =  0.04,  Sct  =  1,  S  =  y/SijSij  and  Ax,  is  the  characteristic  length  of  the  filter.  This 
model  considers  the  anisotropic  part  of  the  SGS  stress  tensor  aXj  =  TL,(uixUj)  —  2/3  k  Stj. 
The  isotropic  components  are  absorbed  in  the  pressure  field. 

For  comparison,  the  DNS  data  are  transposed  from  the  original  high  resolution  193s  points 
to  the  coarse  333  points.  In  the  comparisons,  we  also  consider  the  “resolved”  and  the  “total” 
components  of  the  Reynolds  averaged  moments.  The  former  are  denoted  by  R(a,  b )  with 
R(a,b)  =  ^(o)  —  (a)^  ^(6)  —  (6)^;  and  the  latter  is  r(o,  6)  with  r(a,6)  =  (a  —  a)  (b  —  b).  In 
DNS,  the  “total”  SGS  components  are  directly  available,  while  in  LES  they  are  approximated 
by  r(a,b)  ~  R(a,b )  +  r(a,b)  [44].  Unless  indicated  otherwise,  the  values  of  the  model 
constants  are  Co  =  2.1,  Ce  =  1,  Cq  —  1;  but  the  effects  of  these  parameters  on  the 
predicted  results  are  assessed. 

Figure  11  show  the  instantaneous  iso-surface  of  the  (4>)  =  0.5  field  t  =  80.  By  this 
time,  the  flow  has  gone  through  pairings  and  exhibits  strong  3D  effects.  This  is  evident  by 
the  formation  of  large  scale  span-wise  rollers  with  presence  of  mushroom  like  structures  in 
streamwise  planes  [45].  Similar  to  previous  results  [16],  the  amount  of  SGS  diffusion  with 
the  Smagorinsky  model  is  significant.  Thus,  the  predicted  results  are  overly  smooth.  The 
Reynolds  averaged  values  of  the  filtered  scalar  field  at  t  =  80  are  shown  in  Fig.  12,  and  the 
temporal  variation  of  the  “scalar  thickness,” 


m  =  \vm= o-9)l  +  |v«*>= o.i) 


(33) 


is  shown  in  Fig.  13.  The  filtered  and  unfiltered  DNS  data  yield  virtually  indistinguishable 
results.  The  dissipative  nature  of  the  Smagorinsky  model  at  initial  times  resulting  in  a 
slow  growth  of  the  layer  is  shown.  All  VSFDF  predictions  compare  well  with  DNS  data  in 
predicting  the  spread  of  the  layer. 
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Several  components  of  the  planar  averaged  values  of  the  second  order  SGS  moments 
are  compared  with  DNS  data  in  Figs.  14-15  for  several  values  of  the  model  constants.  In 
general,  the  VSFDF  results  are  in  better  agreement  with  DNS  data  than  those  predicted 
by  the  Smagorinsky  model.  In  this  regard,  therefore,  the  VSFDF  is  expected  to  be  more 
effective  than  the  Smagorinsky  type  closures  for  LES  of  reacting  flows  since  the  extent  of 
SGS  mixing  is  heavily  influenced  by  these  SGS  moments  [52,  53].  However,  it  is  not  possible 
to  suggest  “optimum”  values  for  the  model  constants,  except  that  at  small  Ce  and  Cq  values, 
the  SGS  energy  is  very  large. 

Several  components  of  the  resolved  second  order  moments  are  presented  in  Figs.  16-17.  As 
expected,  the  performance  of  the  Smagorinsky  model  is  not  very  good  as  it  does  not  predict 
the  spread  and  the  peak  value  accurately.  The  VSFDF  yields  reasonable  predictions  except 
for  small  Cq  values.  However,  the  total  values  of  these  moments  are  fairly  independent  of 
the  model  constants  and  yield  very  good  agreement  with  DNS  data  as  shown  in  Figs.  18-19. 
It  is  also  noted  that  while  the  SGS  moments  and/or  the  resolved  moments  may  be  over- 
and/or  under-estimated  depending  on  the  values  of  the  model  coefficients,  the  total  values 
of  the  moments  are  fairly  independent  of  these  coefficients,  at  least  in  the  range  of  values 
as  considered.  But  low  values  of  Cq,  Ce  are  not  recommended  as  they  would  result  into  too 
much  SGS  energy  in  comparison  to  the  resolved  energy. 

VI.  SUMMARY  AND  CONCLUDING  REMARKS 

The  filtered  density  function  (FDF)  methodology  has  proven  effective  for  LES  of  turbulent 
reactive  flows.  In  previous  investigations,  either  the  marginal  FDF  of  the  scalar,  or  that  of 
the  velocity  were  considered.  The  objective  of  present  work  is  to  develop  the  joint  velocity- 
scalar  FDF  methodology.  For  this  purpose,  the  exact  transport  equation  governing  the 
evolution  of  VSFDF  is  derived.  It  is  shown  that  effects  of  the  SGS  convection  and  chemical 
reaction  appear  in  a  closed  form.  The  unclosed  terms  are  modeled  in  a  fashion  similar 
to  those  typically  followed  in  PDF  methods.  The  modeled  VSFDF  transport  equation 
is  solved  numerically  via  a  Lagrangian  Monte  Carlo  (MC)  scheme  via  consideration  of  a 
system  of  equivalent  stochastic  differential  equations  (SDEs).  These  SDEs  are  discretized 
via  the  Euler-Maruyamma  approximation.  The  consistency  of  the  VSFDF  method  and  the 
convergence  of  its  MC  solutions  are  assessed  in  LES  of  a  two-dimensional  (2D)  temporally 
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developing  mixing  layer.  This  assessment  is  done  by  comparing  the  results  obtained  by  the 
MC  procedure  with  those  of  the  finite-difference  scheme  (LES-FD)  for  the  solution  of  the 
transport  equations  of  the  first  two  moments  of  VSFDF.  By  including  the  third  moments 
from  the  VSFDF  into  the  LES-FD,  the  consistency  and  convergence  of  the  MC  solution  are 
demonstrated  by  good  agreements  of  the  first  two  SGS  moments  with  those  obtained  by 
LES-FD. 

The  VSFDF  predictions  are  compared  with  LES  results  with  the  Smagorinsky  [18]  SGS 
model.  All  of  these  results  are  also  compared  with  direct  numerical  simulation  (DNS)  data 
of  a  three-dimensional,  temporally  developing  mixing  layer.  It  is  shown  that  the  VSFDF 
performs  well  in  predicting  some  of  the  phenomena  pertaining  to  the  SGS  transport.  Most  of 
the  overall  flow  features,  including  the  mean  field,  the  resolved  and  total  stresses  as  predicted 
by  VSFDF  are  in  good  agreement  with  DNS  data.  It  may  be  possible  to  improve  the 
predictive  capabilities  of  the  VSFDF  by  two  ways:  (1)  development  of  a  dynamic  procedure 
to  determine  the  model  coefficients,  and/or  (2)  implementation  of  higher  order  closures 
for  the  generalized  Langevin  model  parameter  Gij  [50].  Future  work  is  recommended  for 
development  and  application  of  VSFDF  for  LES  of  chemically  reactive  turbulent  flows. 
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Captions 


Table  1.  Attributes  of  the  computational  methods. 

Figure  1.  Concept  of  ensemble  averaging.  Shown  are  three  different  ensemble  domains: 
1(A#  =  A/2,  Ne  ~  10),  2(A e  —  A,  Ne  «  40),  3(A e  =  2A,  Ne  ~  160).  Black  squares  de¬ 
note  the  finite-difference  grid  points,  and  the  circles  denote  the  MC  particles. 

Figure  2.  Cross-stream  variation  of  the  Reynolds-averaged  values  of  { p )  at  t=34.3:  (a) 
Ne  =  40,  (b)  A E  =  A/2. 

Figure  3.  Temporal  evolution  of  the  scalar  (with  superimposed  vorticity  iso-lines)  (left) 
and  the  vorticity  (right)  fields  for  LES-FD,  with  A e  =  A/2  and  Ne  =  40  at  several  times. 

Figure  4.  Temporal  evolution  of  the  scalar  (with  superimposed  vorticity  iso-lines)  (left) 
and  the  vorticity  (right)  fields  for  VSFDF  with  A#  —  A/2  and  Nb  —  40  at  several  times. 

Figure  5.  Statistical  variability  of  LES-FD  and  VSFDF-C  simulations  with  Ne  =  40  for 
Reynolds-averaged  values  of  r(u,(f>)  at  t=34.4.  Thick  lines:  LES-FD,  thin  lines:  VSFDF-C. 

Figure  6.  Statistical  variability  of  LES-FD  and  VSFDF-C  simulations  with  Ne  =  40  for 
Reynolds-averaged  values  of  t(v,  <j) )  at  t=34.4.  Thick  lines:  LES-FD,  thin  lines:  VSFDF-C. 

Figure  7.  Cross-stream  variations  of  the  Reynolds-averaged  values  of  r(u,(f>)  (a)  A e  = 
A/2,  (b)  A e  —  A,  (c)  Aj 5  =  2A. 

Figure  8.  Cross-stream  variations  of  the  Reynolds-averaged  values  of  t(v,4> )  (a)  A#  = 
A/2,  (b)  A E  =  A,  (c)  A E  =  2A. 

Figure  9.  Cross-stream  variations  of  the  Reynolds-averaged  values  of  r(u,  4>)  (a)  Ne  =  20, 
(b)  Ne  =  40,  (c)  Ne  =  80. 

Figure  10.  Cross-stream  variations  of  the  Reynolds-averaged  values  of  t(v,  <f>)  (a)  Ne  = 
20,  (b)  Ne  =  40,  (c)  NE  =  80. 

Figure  11.  Contours  surface  of  the  ( <f> )  field  in  the  3D  mixing  layer  at  t  =  80  as  obtained 
by:  (a)  DNS,  (b)  Smagorinsky,  (c)  VSFDF. 

Figure  12.  Cross-stream  variations  of  the  Reynolds  averaged  values  of  the  filtered  scalar 
field  at  t  =  80. 

Figure  13.  Temporal  variations  of  the  scalar  thickness. 

Figure  14.  Cross  stream  variations  of  some  of  the  components  of  r  at  t  =  60. 

Figure  15.  Cross  stream  variations  of  some  of  the  components  of  r  at  t  =  80. 

Figure  16.  Cross-stream  variations  of  some  of  the  components  of  R  at  t  =  60. 
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Figure  17.  Cross-stream  variations  of  some  of  the  components  of  R  at  t  —  80. 
Figure  18.  Cross  stream  variations  of  r  at  t  =  60. 

Figure  19.  Cross  stream  variations  of  F  at  t  =  80. 
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(a)  t=13.98  (b)  t=23.94  (c)  t=33.90 


13.98  (b)  t=23.94  (c)  t=33.90 
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Appendix  III 


This  Appendix  is  published  as  Ref.  [21]. 
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Abstract.  An  overview  is  presented  of  recent  developments  and  contribu¬ 
tions  in  large  eddy  simulation  (LES)  of  turbulent  reactive  flows.  The  foun¬ 
dation  of  some  of  the  recently  proposed  subgrid  scale  (SGS)  closures  for 
such  simulations  is  presented,  along  with  a  discussion  of  their  capabilities 
and  limitations.  The  scope  of  the  review  is  limited  to  physical  modeling. 
In  doing  so,  only  issues  pertaining  to  additional  complexities  caused  by 
chemical  reactions  are  discussed.  That  is,  the  challenges  associated  with 
“general”  LES  of  non-reactive  flows  are  not  considered,  even  though  all 
of  these  challenges  are  indeed  present  (and  in  most  cases  are  a  lot  more 
complex)  in  reactive  flows.  It  is  recognized  that  numerical  algorithms  and 
computational  procedures  play  a  significant  role  in  (any)  LES.  However, 
this  review  does  not  deal  with  these  issues  except  for  cases  wherein  the 
actual  numerical-computational  methodology  is  directly  coupled  with  the 
procedure  by  which  LES  is  conducted.  The  SGS  closure  based  on  the  re¬ 
cently  developed  “filtered  density  function”  (FDF)  method  is  described  in 
a  greater  detail.  This  is  due  to  more  familiarity  of  this  reviewer  with  this 
closure;  it  does  not  imply  that  other  closures  are  less  effective. 


1.  Introduction 

In  the  late  1980?s,  I  was  preparing  a  review  article  on  large  scale  numerical 
simulations  of  turbulent  reactive  flows.  The  intent  was  to  provide  a  sur¬ 
vey  of  the  contributions  made  to  both  direct  numerical  simulation  (DNS) 
and  large  eddy  simulation  (LES).  However,  when  that  article  was  finally 
published  (Givi,  1989),  its  content  was  heavily  biased  towards  DNS.  This 
was  not  intentional,  it  just  reflected  the  state  of  progress  on  LES  of  turbu- 
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lent  combustion  at  that  time.  But  with  all  of  the  enthusiasm  for  DNS  in 
the  combustion  community,  the  limitations  of  such  simulations  were  well 
recognized  (even  with  the  most  optimistic  predictions  of  growth  in  super¬ 
computer  technology).  It  was  also  clear  that  the  future  of  large  scale  sim¬ 
ulations  of  practical  turbulent  reacting  flows  would  heavily  depend  on  the 
development  of  LES.  Therefore,  it  was  quite  easy  to  predict  that  LES  would 
receive  significant  attention  in  computational  predictions  of  turbulent  re¬ 
acting  flows  in  the  1990’s  and  into  the  next  (present)  century. 

Now,  at  the  time  of  writing  this  article  (Summer  2001),  while  struggling 
to  meet  the  deadline  for  its  submission(I)  I  am  not  surprised  by  the  extent 
of  the  contributions  in  developing  subgrid  scale  (SGS)  models  or  by  the 
magnificent  work  on  LES  of  a  variety  of  turbulent  reacting  flow  systems. 
In  fact,  I  admit  that  the  rate  of  these  developments  has  been  a  lot  faster 
than  my  capability  to  absorb,  or  in  some  cases  even  follow,  the  details  of 
the  proposed  methodologies.  In  addition,  the  page-limit  restrictions  under 
which  this  is  being  prepared,  preclude  describing  the  details  of  the  wide 
variety  of  currently  available  closures;  similarly,  citation  of  the  relevant 
references  cannot  even  be  done  exhaustively.  Fortunately,  many  aspects  of 
SGS  closures  and  LES  of  reacting  turbulence  have  recently  been  discussed  in 
several  excellent  tutorial  and  review  articles  (Cook  and  Riley,  1998a;  Candel 
et  a/.,  1999;  Bilger,  2000;  Branley  and  Jones,  2000;  Menon,  2000;  Peters, 
2000;  Pope,  2000;  Luo,  2001;  Poinsot  and  Veynante,  2001).  Therefore,  in 
the  present  review  I  concentrate  on  some  of  the  major  issues  related  to  my 
area  of  research  within  this  field. 

2.  Starting  Equations 

Large  eddy  simulation  involves  the  use  of  the  spatial  filtering  operation 
(Sagaut,  2001) 


/+oo 

<2(x',  <)£(*'>  x)dx',  (l) 

-oo 

where  Q  denotes  the  filter  function  of  width  Aq,  and  (£J(x,£))^  represents 
the  filtered  value  of  the  transport  variable  Q(x,  t ).  In  variable  density  flows 
it  is  convenient  to  consider  the  Favre  filtered  quantity, 

(<2(x,  t))L  =(pQ)e/(p}e.  (2) 

We  consider  spatially  &  temporally  invariant  and  localized  filter  functions, 
t/(x',x)  =  G(x'-x)  with  the  properties  G(x)  =  G(-x),  and  G(x)dx  = 
1.  Moreover,  we  only  consider  “positive”  filter  functions  for  which  all  the 
moments  xmG(x)dx  exist  for  m  >  0. 
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To  set  the  framework,  we  consider  the  transport  equations  of  chemically 
reacting  flows.  To  isolate  the  effects  of  chemical  reaction  in  the  simplest 
way,  we  consider  single-phase  (gaseous)  combustion  in  a  low  Mach  number 
flow  with  negligible  radiative  heat  transfer  and  viscous  dissipation.  We  also 
assume  that  Newton’s  law  of  viscosity,  Fourier’s  law  of  heat  conduction  and 
Fick’s  law  of  mass  diffusion  are  applicable.  Therefore,  the  primary  transport 
variables  are  the  density  p,  the  velocity  vector  i  =  1,2,3  along  the 
X{  direction,  the  pressure  p.  the  species’  mass  fractions  Ya,  and  the  total 
specific  enthalpy  h.  All  of  the  mass  fractions  and  the  enthalpy  are  grouped 
into  the  scalar  array  <£(x,t)  =  [4>i,4>2,  ■■■4>o]  =  [ii,  Y2, . .  .,Yns,Ii]  of  size 
a  —  Ns  +  1  where  Ns  denotes  the  total  number  of  species.  Application  of 
the  filtering  operation  to  the  equations  of  continuity,  momentum,  enthalpy 
(energy)  and  species  mass  fraction  equations  gives 

d(p)t  .  d(p)t{v,i)L  _ 

dt  +  dxi  -u’  w 

mt(»j)L  ,  d(p)e(ui)L(uj)L  d(p)t  d(Tjj)t  dTjj 

dt  dxi  dxj  dxi  dxi  ’  '  ' 

— m — + — w< —  ~  ~  + ipS°}‘'  <5) 

where  t  represents  time,  and  the  filtered  reaction  source  terms  are  denoted 
by  (pSa}i  —  The  viscous  stress  tensor  and  the  mass/heat  fluxes 

are  denoted  by  and  Jf ,  respectively.  At  low  Mach  numbers  and  heat 
release  rates,  by  neglecting  the  viscous  dissipation  and  thermal  radiation 
the  source  terms  in  the  enthalpy  equation  can  be  assumed  to  be  negligible. 
Thus,  SQ  =  Sa(<f> ).  The  terms  Tij  =  (p)t((uiUj)L  -  ( Ui)L(uj)L )  and  Mf  = 
{p)i{{ui<t>a)L  ~  {^i) L(<t>a) l)  denote  the  SGS  stress  and  the  SGS  mass  flux, 
respectively.  Equations  (3)-(5)  are  coupled  through  the  equation  of  state. 

3.  Closure  Methodologies 

For  non-reacting  flows  the  SGS  closure  is  associated  with  Tij  and  M/* 
(Canuto,  1994;  Ciofalo,  1994;  Lesieur  and  Metais,  1996).  In  reacting  flows, 
an  additional  model  is  required  for  the  filtered  reaction  rate  (S^l-  This 
modeling  is  the  subject  of  primary  concern  in  this  review. 

One  of  the  first  contributions  in  LES  of  reactive  flows,  similar  to  that  in 
LES  of  non-reacting  flows,  was  made  in  atmospheric  sciences  (Schumann, 
1989).  In  this  work,  the  effects  of  SGS  scalar  fluctuations  (as  appear  in  the 
chemical  source  term)  are  assumed  negligible,  i.e.  (Sa{<f>))L  ^  tow. 
This  assumption  is  compatible  with  that  made  in  some  of  the  more  recent 
contributions  (Boris  et  ai,  1992;  Fureby  and  Grinstein,  1999),  in  which  it 
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is  argued  that  all  of  the  essential  SGS  contributions  are  included  in  the 
numerical  discretization  procedure. 

Modeling  of  the  scalar  fluctuations  has  been  the  subject  of  broad  in¬ 
vestigations  in  Reynolds  averaged  simulations  (RAS)  for  over  five  decades, 
resulting  in  a  variety  of  closure  strategies  (Libby  and  Williams,  1980;  Libby 
and  Williams,  1994).  Within  the  past  10  years  or  so,  almost  all  of  these 
closures  have  been  considered  for  LES.  Examples:  the  eddy- break  up  mod¬ 
els  (Fureby  and  Lofstrom,  1994;  Candel  et  al.,  1999),  moment  methods 
(Frankel  et  al.,  1993),  the  flamelet  concept  (Cook  et  al.,  1997;  Cook  and 
Riley,  1998b;  De  Bruyn  Kops  et  al.,  1998;  DesJardin  and  Frankel,  1998; 
DesJardin  and  Frankel,  1999;  Pitsch  and  Steiner,  2000;  Ladeinde  et  al., 
2001),  the  linear  eddy  model  (LEM)  (McMurtry  et  al.,  1992;  Menon  and 
Calhoon,  1996;  Kim  et  al.,  1999;  Menon,  2000),  the  conditional  moment 
method  (CMM)  (Bushe  and  Steiner,  1999;  Steiner  and  Bushe,  2001),  and 
many  others  (Sykes  et  al.,  1992;  Galperin  and  Orszag,  1993;  Smith  and 
Menon,  1996;  Im  et  al.,  1997;  McGrattan  et  al,  1998;  Thibaut  and  Candel, 
1998;  Battaglia  et  al.,  2000;  Collin  et  al.,  2000).  In  addition,  several  of  the 
closures  previously  developed  for  LES  of  non-reacting  flows,  have  been  ex¬ 
tended  for  use  in  reacting  flow  simulations  (DesJardin  and  Frankel,  1998; 
Jaberi  and  James,  1998). 

The  probability  density  function  (PDF)  methods  have  proven  particu¬ 
larly  useful  in  RAS  (O’Brien,  1980;  Pope,  1985;  Dopazo,  1994;  Fox,  1996; 
Pope,  2000).  The  systematic  approach  for  determining  the  PDF  is  by  means 
of  solving  its  transport  equation.  An  alternative  approach  is  based  on 
assumed  methods  in  which  the  shape  of  the  PDF  is  specified  a  priori. 
This  has  been  pursued  in  several  studies  in  most  of  which  it  is  assumed 
that  the  thermo-chemical  variables  depend  only  on  the  mixture  fraction, 
e.g.  infinitely  fast  reaction,  equilibrium  chemistry.  Therefore,  the  PDF  is 
univariate  (Madnia  and  Givi,  1993;  Cook  and  Riley,  1994;  Reveillon  and 
Vervisch,  1996;  Branley  and  Jones,  1997;  Jimenez  et  al.,  1997;  Mathey  and 
Chollet,  1997;  DesJardin  and  Frankel,  1998;  DesJardin  and  Frankel,  1999; 
Forkel  and  Janicka,  2000;  Kempf  et  al.,  2000).  For  LES  of  non-equilibrium 
reactive  flows,  it  is  necessary  to  assume  the  joint  PDF  of  multi-scalars 
(Frankel  et  al.,  1993).  Consistent  with  popular  methods  of  generating  uni¬ 
variate  (Leemis,  1986)  and  multivariate  (Johnson  and  Kotz,  1972)  distribu¬ 
tions,  all  of  the  assumed  SGS  scalar  PDFs  in  the  contributions  cited  above 
are  based  on  the  first  and  the  second  order  moments.  The  PDFs  generated 
in  this  way  offer  sufficient  flexibility  and  are  affordable  for  large  scale  sim¬ 
ulations.  However,  it  is  now  well  understood  that  the  “true”  PDF  strongly 
depends  on  the  actual  physics  of  mixing  in  a  given  flow  condition  (Jaberi 
et  al.,  1996).  Therefore,  there  is  a  need  to  determine  such  PDFs  in  a  more 
systematic  manner. 
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The  “filtered  density  function”  (FDF)  methodology  introduced  by  Pope 
(1990)  provides  the  framework  for  fundamental  developments  of  the  PDF 
based  SGS  closures.  This  method  provides  a  means  of  determining  the  PDF 
from  its  own  transport  equation.  For  the  scalars’  array  <p(x,t )  the  FDF, 
denoted  by  Pl,  is  defined  as  (Pope,  1990) 

/+0° 

C  bl>,  <A(x',  *)]  G(xr  -  x)dx',  (6) 
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C  [V>,  <Kx,  <)]  =  n  ~  &»(x>  *)]>  (7) 

a= 1 

where  6  denotes  the  delta  function  and  ip  denotes  the  composition  domain 
of  the  scalar  array.  The  term  ([(p  —  ip(x,  t)]  is  the  “fine-grained”  density 
(Lundgren,  1967;  O’Brien,  1980;  Pope,  1985;  Dopazo,  1994).  In  variable 
density  flows,  it  is  convenient  to  consider  the  “filtered  mass  density  func¬ 
tion”  (FMDF),  denoted  by  Fl,  as 


/+oo 

Kx',  OC  Ns  <£(x',  t)\  G(x'  -  x)dx!.  (8) 
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The  integral  property  of  the  FDF  and  FMDF  is  such  that 

/•foo  r+oo 

PL(ip;x,t)dip  =  1,  /  FL(ip;x,t)dip  =  (p(x,t))e.  (9) 

-oo  J — oo 

For  further  discussions,  it  is  useful  to  define  the  mass  weighted  conditional 
filtered  mean  of  the  variable  Q(x,t), 


_  f-< »  t)Q(x',  t)(  [ip,  <p(x',  t)]  G(x!  -  x)dx' 


(Q(x,t)\ip)t  = 


FL(ip;x,t) 


(10) 


Therefore,  when  Q  can  be  completely  described  by  the  compositional  vari¬ 
able,  i.e.  Q(x,t)  =  Q(0(x,t)),  we  have  =  Q{^).  Also, 

/+oo 

(<9(x,  t)\ip)eFL(ip ;  x,  t)dip  =  (p(x,  t))t{Q(x,  t))L.  (11) 

-oo 


The  transport  equation  for  Fi{ip]  x,  /)  is  obtained  by  multiplying  the  trans¬ 
port  equation  for  the  fine  grained  density  by  the  filter  function  G(x!  —  x) 
and  integrating  over  x'  space  (Gao  and  O’Brien,  1993;  Colucci  et  al.,  1998; 
Reveillon  and  Vervisch,  1998;  Jaberi  et  al.,  1999;  Jaberi,  1999;  Zhou  and 
Pereira,  2000;  Tong,  2001), 
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dFi,(ifr;x,t)  x,  <)] 

dt  dxi 


d[Sa(VQFjr,(V>;x>*)] 

dlpa 


[<a 


The  first  term  on  the  RHS  is  due  to  chemical  reaction  and  is  in  a  closed 
form.  This  demonstrates  the  primary  advantage  of  the  FDF  methodology. 
However,  the  SGS  convection  (the  second  term  on  the  LHS)  and  SGS  mix¬ 
ing  (the  second  term  on  the  RHS)  must  be  modeled.  One  of  the  most  chal¬ 
lenging  issues  in  FDF  is  associated  with  closure  of  the  mixing  term.  This 
has  been  the  subject  of  broad  investigations  in  PDF  modeling  (Pope,  1985; 
Pope,  2000).  In  Eq.  (12)  the  effects  of  mixing  are  displayed  through  the 
“conditional  expected  diffusion”  of  the  scalars,  but  can  also  be  represented 
in  the  form  of  the  “conditional  expected  dissipation”  (O’Brien,  1980;  Pope, 
1985).  The  closure  for  this  can  be  via  any  of  the  ones  currently  in  use  in 
PDF  methods  (Pope,  2000).  In  the  absence  of  a  clearly  superior  model,  the 
linear  mean  square  estimation  (LMSE)  model  (O’Brien,  1980)  has  been 
used  in  almost  all  of  previous  LES  based  on  FDF  (Colucci  et  al.,  1998; 
Jaberi  et  al.,  1999;  Garrick  et  al.,  1999;  James  and  Jaberi,  2000;  Zhou  and 
Pereira,  2000).  With  Jf  =  — 7§^S  this  model  is 


d  \d(FL/(p)t) 


T  n(V,a  {fic^L^Fl,]  ,(13) 


where  flm(x,  t)  is  the  “frequency”  of  mixing  within  the  subgrid  and  must 
be  modeled.  The  convective  term  can  be  modeled  as 


(uiWeFL  =  (u^lFl  -  7 1 


d(FL/(p)t ) 


where  7 <  is  the  SGS  diffusion  coefficient  and  must  be  specified.  Equation 
(14)  is  in  accord  with  that  often  used  in  conventional  LES  (Moin  et  al., 
1991;  Canuto,  1994;  Ciofalo,  1994;  Lesieur  and  Metais,  1996).  With  this 
formulation,  obviously  the  resolved  hydrodynamic  field  must  be  determined 
by  other  means.  This  problem  can  be  circumvented  by  considering  the  joint 
velocity-scalar  FMDF, 
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3  <7 

f  [v,  u(x,  t),  I)]  =  JJ  6[vk  -  uk(x, <)]  (16) 

k= 1  ot= 1 

where  v  denotes  the  composition  domain  of  the  random  velocity  vector,  and 
£[v,  u(x,t),V>,  <£(x,f)]  is  the  fine-grained  velocity-scalar  density.  Most  re¬ 
cent  work  in  this  regard  consider  the  transport  of  the  velocity  FDF  (VFDF) 
(Gicquel,  2001)  and  the  joint  velocity-scalar  FDF  (VSFDF)  (Drozda,  2001). 
The  operational  procedure  is  similar  to  that  developed  previously  for  PDF 
methods  (Pope,  1985;  Pope,  1994;  Pope,  2000). 

The  closure  problems  as  noted  above  are  not  particular  to  the  FDF; 
all  of  the  other  schemes  require  similar  modelings.  For  example,  in  the 
limit  of  equilibrium  chemistry  all  of  the  statistics  of  the  reacting  fields 
are  related  to  those  of  the  mixture  fraction.  The  FMDF  of  the  mixture 
fraction  can  be  obtained  from  the  solution  of  Eq.  (12)  with  S  =  0.  So, 
there  is  still  a  need  for  modeling  of  the  mixing  term.  Even  in  cases  where 
the  FMDF  is  assumed,  its  distribution  is  parameterized  with  the  low  order 
moments  of  the  mixture  fraction.  As  indicated  above,  the  first  two  moments 
are  typically  used  for  this  parameterization.  Therefore,  there  is  a  need  for 
closure  of  the  “total  SGS  dissipation”  as  appears  in  the  second  moment 
(SGS  variance)  equation.  Several  means  of  dealing  with  this  closure  problem 
are  available  (Girimaji  and  Zhou,  1996;  Pierce  and  Moin,  1998;  Cook  and 
Bushe,  1999;  Jimenez  et  al.,  2001;  De  Bruyn  Kops  and  Riley,  2001). 

The  above  problem  is  a  bit  more  complex  when  the  SGS  chemical  reac¬ 
tion  is  assumed  to  be  in  the  “flamelet”  regime  (Peters,  2000).  In  this  case, 
even  with  the  one-dimensional  flamelet  model,  the  thermo- chemical  vari¬ 
ables  are  parameterized  by  the  mixture  fraction  and  its  rate  of  dissipation 
(Cook  et  al.,  1997;  Cook  and  Riley,  1998b;  De  Bruyn  Kops  et  al.,  1998; 
DesJardin  and  Frankel,  1998;  Cook  and  Riley,  1998b).  Therefore,  there  is  a 
need  for  a  priori  specification  of  the  joint  FDF  of  the  mixture  fraction  and 
its  dissipation.  A  review  of  different  methods  of  dealing  with  this  issue  is 
available  (Cook  and  Riley,  1998a).  Equation  (12)  with  5  =  0  indicates  that 
there  is  a  dependence  between  the  FDF  of  the  mixture  fractions  and  the 
conditional  expected  diffusion  (and  the  conditional  expected  dissipation). 
This  dependency  is  not  considered  in  most  previous  contributions,  but  is 
the  subject  of  current  investigations  (DesJardin  et  al.,  2001). 

Modeling  of  the  conditional  expected  dissipation  is  also  required  in  the 
conditional  moment  method  (Bushe  and  Steiner,  1999;  Steiner  and  Bushe, 
2001).  This  issue  has  been  recognized  at  the  early  stages  of  developments  of 
CMM  in  RAS  (Bilger,  2000).  With  this  model,  the  conditional  filtered  mean 
values  of  the  thermo-chemical  variables  (LHS  of  Eq.  (10))  are  obtained 
by  their  modeled  transport  equation.  This  is  obviously  computationally 
less  demanding  that  solving  the  FDF  transport  equation.  But  in  order  to 
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determine  the  actual  filtered  quantities,  the  distribution  of  the  mixture 
fraction  FDF  must  be  specified. 

An  important  issue  in  regard  to  FDF  is  associated  with  the  numeri¬ 
cal  solution  of  its  transport  equation.  The  Lagrangian  Monte  Carlo  scheme 
(Pope,  1985)  has  proven  particularly  useful  for  this  purpose.  In  this  scheme, 
the  FDF  is  represented  via  an  ensemble  of  computational  elements  or  par¬ 
ticles.  Transport  of  these  particles  and  the  change  in  their  properties  are 
modeled  by  a  set  of  stochastic  differential  equations  (SDEs)  (Soong,  1973). 
The  diffusion  process  (Gardiner,  1990)  has  proven  effective  for  this  pur¬ 
pose.  The  coefficients  in  the  Langevin  equation  governing  this  process  are 
set  in  such  a  way  that  the  resulting  Fokker- Planck  equation  (Risken,  1989) 
is  equivalent  to  the  FDF  transport  equation.  Therefore,  the  Monte  Carlo 
solution  of  the  SDEs  represent  the  solution  of  the  FDF  in  the  probabilistic 
sense.  This  procedure  has  proven  successful  for  simulating  PDF  in  a  variety 
of  systems  (Grigoriu,  1995).  However,  one  must  be  careful  in  performing 
stochastic  simulations  in  conjunction  with  modern  CFD  solvers.  Many  of 
the  advanced  discretization  routines  developed  for  solving  deterministic  dif¬ 
ferential  equations  may  not  be  applicable,  or  may  have  to  be  significantly 
modified  to  be  suitable  for  solving  SDEs  (Kloeden  and  Platen,  1995). 

Implementation  of  LEM  is  also  based  on  stochastic  representation  of 
the  flow.  In  its  original  development  in  RAS  (Kerstein,  1988),  the  processes 
of  molecular  diffusion,  chemical  reaction  and  turbulent  convection  are  con¬ 
sidered  separately.  This  is  achieved  by  a  reduced  one-dimensional  (linear) 
description  of  the  scalar  field,  which  makes  it  possible  to  resolve  the  flow 
scales  even  for  flows  with  relatively  high  Reynolds,  Schmidt  and  Damkohler 
numbers.  The  interpretation  of  the  one-dimensional  domain  is  dependent 
on  the  particular  flow  under  consideration.  In  this  way,  the  processes  of 
molecular  diffusion  and  chemical  reaction  are  taken  into  account  exactly, 
but  the  effects  of  convection  are  modeled.  This  is  achieved  by  “random 
rearrangement”  (or  stirring)  events  in  such  a  way  that  the  displacements 
of  fluid  elements  result  in  a  diffusivity  equal  to  the  “turbulent  diffusivity.” 
For  LES,  this  procedure  is  followed  within  each  of  the  computational  cells, 
and  stirring  is  performed  to  yield  the  desired  SGS  diffusivity.  Menon  and 
colleagues  have  made  extensive  use  of  LEM  for  LES  of  a  wide  variety  of 
reacting  flows.  A  recent  review  is  available  (Menon,  2000). 
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